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We give a detailed analysis of the origin of spurious divergences and finite steps that have been 
recently identified in particle-number restoration calculations within the nuclear energy density func- 
tional framework. We isolate two distinct levels of spurious contributions to the energy. The first 
one is encoded in the definition of the basic energy density functional itself whereas the second one 
relates to the canonical procedure followed to extend the use of the energy density functional to 
multi-reference calculations. The first level of spuriosity relates to the long-known self-interaction 
problem and to the newly discussed self-pairing interaction process which might appear when de- 
scribing paired systems with energy functional methods using auxiliary reference states of Bogoli- 
ubov or BCS type. A minimal correction to the second level of spuriosity to the multi-reference 
nuclear energy density functional proposed in [D. Lacroix, T. Duguet, M. Bender, arXiv:0809.204l] 
is shown to remove completely the anomalies encountered in particle-number restored calculations. 
In particular, it restores sum-rules over (positive) particle numbers that are to be fulfilled by the 
particle-number-restored formalism. The correction is found to be on the order of several hundreds 
of keVs up to about 1 MeV in realistic calculations, which is small compared to the total binding 
energy, but often accounts for a substantial percentage of the energy gain from particle-number 
restoration and is on the same energy scale as the excitations one addresses with multi-reference 
energy density functional methods. 

PACS numbers: 21.10.Rc, 21.60.Ev, 71.15.Mb 



I. INTRODUCTION 

Methods based on the use of energy density function- 
als (EDF) [l| currently provide the only set of theoretical 
tools that can be applied to all nuclei but the lightest ones 
in a systematic manner irrespective of their mass and 
isospin. Nuclear EDF methods coexist on two distinct 
levels. On the first level, that is traditionally and inap- 
propriately called "self-consistent mean- field theory" or 
Hartree-Fock (HF) or Hartree-Fock-Bogoliubov (HFB) , a 
single product state provides the normal and anomalous 
density matrices the energy is a functional of. Wc will 
call this type of method a singlc-rcfcrcncc (SR) EDF ap- 
proach. On the second level, traditionally and inappro- 
priately called " beyond- mean- field methods", i.e. sym- 
metry restoration and configuration mixing in the spirit 
of the Generator Coordinate Method (GCM), the set 
of transition density matrices defined from an appro- 
priate ensemble of product states enter the EDF. We 
will call such a method a multi-reference (MR) EDF ap- 
proach. Although SR EDF calculations have many simi- 
larities with Density Functional Theory (DFT) which is 
widely used in atomic, molecular and condensed matter 
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physics 0, H, H, H, [1, 0, Hi J they also present key differ- 
ences, which prohibit the straightforward mapping of the 
concepts of electronic DFT to the nuclear case 0, [HI ■ 

The reference state entering a SR EDF calculation usu- 
ally breaks several symmetries of the exact eigenstates of 
the nuclear Hamiltonian. This is done on purpose, as 
it allows one to incorporate so-called static correlations 
associated with collective modes [13, [13, [13, [13 mod- 
erate computational cost. One of the most important 
categories of correlations which can be grasped this way 
are those associated with the formation of neutron and 
proton Cooper pairs in the medium. 

In a SR EDF approach, pairing correlations are incor- 
porated by making the energy a functional of the anoma- 
lous density matrix in addition to the normal one. This 
amounts to using an independent quasi-particle state 
(which will be called a quasi-particle vacuum in what fol- 
lows) of BCS or Bogoliubov type as a reference state in- 
stead of a Slater determinant. The price to pay is break- 
ing the U{1) symmetry in gauge space that is a feature of 
eigenstates of the particle-number operator. As a result 
the SR state is spread in particle-number space, and one 
cannot associate the computed energy, even implicitly, 
to a state belonging to a specific irreducible represen- 
tation of U{1). In condensed matter ph ysic s, for which 
the BCS method was originally designed this is usu- 
ally not much of a problem. Nuclei, however, are small 
finite quantum many-body systems for which two prob- 
lems arise in this context: (i) the SR approach docs not 
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grasp the so-called dynamical pairing correlations associ- 
ated with the fluctuations of both the magnitude and the 
phase of the order parameter of the broken t/(l) symme- 
try. Correlations associated with this zero-energy mode 
may affect any observable that probes the occupation 
of levels around the Fermi surface in a significant way; 
(ii) when the density of single-particle levels around the 
Fermi energy is below a critical value, pairing correla- 
tions are entirely dynamical and cannot be described by 
the SR method, in most cases in contradiction with ex- 
periment. 

All of these limitations can be overcome by performing 
multi-reference EDF calculations. Those allow in partic- 
ular the restoration of particle number [TtI. ITsI. fiol. f20l. l2l| . 
It has been noticed for some time, however, that particle- 
number restored energies might exhibit divergences 
H^, IH] and finite steps [13, [23 whenever a single-particle 
level crosses the Fermi energy as a function of a collective 
coordinate. This problem is particular to energy density 
functionals, but absent in approaches based on the use 
of a genuine Hamiltonian and a correlated wave func- 
tion. As pointed out by Anguiano et al. in Ref. p^ . 
some of the common assumptions and approximations 
made in the construction of nuclear EDFs unavoidably 
lead to such anomalies, and these authors, as done ear- 
lier in Refs. [H, US] in a different context, advocate to 
use strict antisymmetric two-body vertices and to keep 
all exchange terms when computing the energy. However, 
and contrary to what is stated in Ref. [isl . , using an- 
tisymmetric but density-dependent two-body vertices is 
not free from pathologies, even when the divergence in- 
troduced by the density-dependent terms is intcgrable. 
There is an additional problem that arises particularly 
when such a dependence is taken under the form of a 
non-integer power of the density (matrix) psi . [27| . 

Practitioners of EDF methods, however, recognize that 
it is desirable to use more general energy functionals. 
For those, Particle-Number Restoration (PNR), and the 
MR formalism in general, still need to be formulated in 
a consistent and unambiguous manner that is free from 
pathologies. As a first step into that direction, a thor- 
ough analysis has been recently given by Dobaczcwski et 
al. regarding the poles and steps contained in a particle- 
number restored energy density functional p5l |. In the 
first of our companion papers |28|, hereafter referred to 
as Paper I, we could connect those pathologies to an un- 
derlying level of spuriosity that is encoded in the SR en- 
ergy functional. The associated spurious terms turn out 
to relate to self-interaction processes well-known in DFT 
for condensed matter [29j . a problem which was actu- 
ally studied beforehand in the nuclear context [sO] but 
was soon forgotten, as well as to spurious self-pairing 
processes, whose notion is introduced in the present pa- 
per. The common source of both pathologies is the 
use of different and non-antisymmetric vertices at dif- 
ferent places in the EDF violating in this way the ex- 
change symmetry of Fermi statistics. The existence of 
spurious self-interaction and self-pairing in the SR en- 



ergy functional is indeed a prerequisite for the appear- 
ance of divergences and steps at the MR level, but it is 
not its origin as such. The pathologies that arc partic- 
ular to the MR level, e.g., particle- number restoration, 
turn out to be caused by an unphysical contribution to 
the weight of the self-interaction and self-pairing con- 
tributions in multi-reference energy kernels. This is an 
unforeseen consequence of the common practice of con- 
structing the multi-reference energy functional kernel by 
replacing the density matrices entering a given SR en- 
ergy functional by transition density matrices [sil . [3^ in 
analogy to the application of the Generalized Wick The- 
orem (GWT) [H, [U within a Hamiltonian- and wave- 
function-based approach. Making reference to a Wick 
theorem in an energy density functional without hav- 
ing a genuine operator to relate to is necessarily out- 
side the scope of that Wick theorem and might pro- 
duce unexpected results. And indeed, using the stan- 
dard [3^] and generalized [H, Wick theorems yields 
different weights to self-interaction and self-pairing con- 
tributions to the MR energy kernel as demonstrated in 
Paper I. Only the GWT-motivated procedure produces 
the poles that are at the origin of the divergences and 
steps, thus introducing a second level of spuriosity. Us- 
ing a Hamiltonian- and wave-function-based approach, 
no problem arises; the vertices at play are either zero or 
recombine in a particular way that cancels out danger- 
ous poles. Our analysis in Paper I was made without 
reference to a particular MR application and aimed at 
the introduction of a proper framework to identify and 
separate both levels of spuriosity within any MR EDF 
calculation. It is the aim of the present paper to apply 
the procedure proposed in Paper I to correct for the un- 
physical weights in the special case of particle-number 
restoration using a particular energy functional the cor- 
rection can be applied to. In a third paper [23|, hereafter 
called Paper III, we analyze in detail in the context of 
PNR the reasons why the pathologies associated with 
more commonly used functionals containing non-integer 
powers of the density (matrix) |25|| ar e very likely to be 
uncorrectable. Together with Ref. [23, Paper III demon- 
strates that the density-dependent two-body forces that 
are advertised by some authors to be free of pathologies 
[isl . [2ll . [26j also have their share of problems when used 
in MR calculations. 

The paper is organized as follows: In Section |TT1 we 
introduce single-reference EDF calculations, paying par- 
ticular attention to resemblances and key differences with 
the HFB method based on the use of a Hamilton opera- 
tor. In Section IIIIl we introduce multi-reference EDF 
calculations appropriate to restoring particle number, 
paying particular attention to resemblances and key dif- 
ferences with the strict particle-number projected HFB 
(PNP-HFB) method based on the use of a Hamilton op- 
erator. In Section HVl we discuss the occurrence of spuri- 
ous self-interaction and self-pairing processes in SR and 
MR calculations. Section |V] analyzes the occurrence of 
spurious self-interaction and self-pairing contributions to 
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the particle-number restored EDF using a complex plane 
analysis and specifies the correction designed in Paper I 
to that particular case. Section IVll applies the correction 
scheme to realistic calculations of finite nuclei. Finally, 
conclusions are drawn in Section FVIII Several appendices 
complement the paper with derivations and formulas use- 
ful for practical applications. 



to choosing a major axis system for quadrupole deformed 
product states. In the case of gauge symmetry, a conve- 
nient orientation is provided by = 0. States at different 
angles are obtained from this state applying the rotation 
operator e"^^ in gauge space 

\%) - e^'^* |<I>o) = e''^* n ("M + < 4) |0) • (2) 



II. SINGLE-REFERENCE EDF APPROACH 

Let us first present the basic elements of the single- 
reference EDF method which will be needed for our dis- 
cussion. The HFB implementation of the single-reference 
EDF approach relies on the use of a quasi-particle vac- 
uum \^ip) as a reference state from which the normal and 
anomalous one-body density matrices entering the energy 
density functional are calculated. In the canonical basis 
{a^,a^} that diagonalizes its one-body normal density 
matrix, the reference state reads 



I*.) = n ' 

fj.>0 



|0>, 



(1) 



where |0) is the particle vacuum. Throughout this pa- 
per we limit ourselves to time-reversal invariant quasi- 
particle vacua |0) with even- number parity and thus only 
discuss explicitly the ground-state of even-even systems. 
In addition, we do not mix protons and neutrons when 
constructing quasi-particle operators. In particular, this 
limits the pairing interaction to particles of the same 
isospin. Identical assumptions are made in most, if not 
all, published work performed using particle-number pro- 
jected energy density functionals so far, and are sufficient 
for the purpose of the present paper. 

The single-particle wave functions associated with the 
pair-conjugated canonical states (/x, p,) is denoted as 0^ 
and A quantum number rj^^ can always be chosen 
to separate the single-particle basis into two halves, the 
"positive" and the "negative" ones, with each partner of 
a given conjugated pair associated to a different half. The 
normalization of |$^) gives -I- \vf^ e^*'^| = 1. We use 
phase conventions where the and Vf^ are real numbers; 



hence, uj^ 



1, which also fixes the global phase of 



\^ip). The angle ip in the remaining phase factor denotes 
the orientation of the state in the U{1) gauge space. 

The exact eigenstates of the nuclear many-body prob- 
lem belong to a specific irreducible representation of the 
U{1) group. By contrast, the product state |$^) behaves 
as a wave packet in gauge space as it mixes states be- 
longing to different irreducible representations. The use 
of such Bogoliubov product states is at the heart of the 
symmetry-breaking description of static pairing correla- 
tions based on a single reference state. In spite of the 
broken symmetry of the product state, all observables 
that are scalars in gauge space still have to be indepen- 
dent on its orientation in gauge space. This allows one 
to choose a convenient angle on the level of SR calcula- 
tions that simplifies the calculations, a procedure similar 



A. Energy in the strict HFB approach 

As a strict HFB approach, we denote the method that 
determines the energetically most favored quasi-particle 
vacuum \^,p) through the minimization of the expecta- 
tion value of a given Hamiltonian H in that product state, 
without any approximations or generalizations. For the 
sake of transparency, the Hamiltonian 



ct Ci 



ijkl 



(3) 



is assumed to be given by the sum of kinetic energy term 
and a two-body interaction. In Eq. ^ {4} defines a 
complete set of single-particle states whereas Vijki de- 
notes antisymmetric matrix elements (or vertices) of the 
two-body interaction in that basis. The discussion be- 
low can be extended without difficulty to a Hamiltonian 
containing three-body or higher-body forces, but this be- 
comes cumbersome and is not necessary for the purpose 
of this paper. 

An important point is that in the context of the strict 
HFB approach, we assume that the vertex Vijki docs not 
depend on density. So-called density-dependent vertices 
of Skyrme and Gogny type are widely used in the liter- 
ature. However, as pointed out in Ref. psj . discussed in 
the present paper and insisted on further in Paper HI, 
any density-dependent effective vertices do provide MR 
energies with (at least) spurious finite contributions, even 
though the vertex is antisymmetric with respect to the 
remaining single-particle degrees of freedom and all as- 
sociated exchange terms are exactly accounted for in the 
MR energy kernels. 

Using the Standard Wick theorem (SWT) [H,!!!,!!!, 
the expectation value of H in the product state can 
be evaluated as 



'^fifiiyy ^fi^fi '^v'^i' 



(4) 



where p'^'^ and k^^ are the normal density matrix and 
anomalous density matrix (pairing tensor) constructed 
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from \^ip), respectively. In the canonical basis of the 
Bogoliubov transformation defining |<i>^), these take the 
simple form 



1/ V p-2»V A _ 



(5) 



(6) 



(7) 



The expectation value given in Eq. ([4]) can be seen as a 
particular functional of p"^"^ , k^"^ and k'^'^*. The sym- 
metries of the Hamiltonian lead of course to a number of 
specific properties of this functional. In particular, since 
the Hamiltonian commutes with the particle-number op- 
erator, one finds that 



E[p^ 



,00 ^00 ^00*1 



(8) 



which underlines that all states that differ only by a ro- 
tation in gauge space are degenerate. In other words, the 
energy functional behaves as a scalar in gauge space as 
expected. 



B. Energy in the SR energy functional approach 

In nuclear physics, strict HFB-type approaches are 
frequently applied in a restricted shell-model space us- 
ing parametrized single-particle energies and an effective 
Hamiltonian as a residual interaction [SS,, j39., .4^. For 
a multitude of reasons outlined in Paper I and refer- 
ences given therein, methods using the full model space 
of occupied particles had to resume so far to the use of 
(phenomenologieal) density-dependent effective interac- 
tions [4l|, which sets the stage for what is nowa- 
days recognized as an approximation to a more gen- 
eral single-reference EDF formalism. This framework 
shares many features with the Density Functional Theory 
(DFT) widely used for description of electronic many- 
body systems 0, H, 0, H, 0, 0, 111 , but also displays key 
differences, which prohibit the straightforward mapping 
of all concepts of electronic DFT to the nuclear case 

BE, [a. 

In the DFT for many-electron systems, constructive 
schemes have been established to design the energy func- 
tional, see for instance Ref. Q and references given 
therein. In nuclear physics, such a procedure that would 
suggest the structure of the functional is still missing, al- 
ready on a qualitative level. The reasons are the complex- 
ity of the nucleon-nucleon interaction on the one hand, 
and that in-medium correlations are never small correc- 
tions on the other hand. In the absence of a construc- 
tive scheme, all widely used nuclear energy functionals 
were set up keeping an underlying two-body and some- 
times three-body interaction as guiding principle, mak- 
ing generalizations suggested by phenomenology and ap- 
proximating or even omitting terms that are small, but 



difficult to evaluate. As a consequence, the structure of 
these functionals resembles that of Eq. (|3]), except that 
the expectation value E [p, k, k*] is replaced by a func- 
tional £ [p, K, K*]. Considering the simple case of a bi- 
linear functional for simplicity and comparison purposes, 
such a functional can be written as 



£ [p, K, 



= £P + £PP + 



' 4 y ^ '^fj.p.fu '^M '^'^ ■ 



(9) 



This might appear as an unsusual way to write standard 
energy functional, but will turn out to be very useful be- 
low. The corresponding explicit expressions for a Skyrme 
energy functional are given in Appendix [X] The crucial 
point for our discussion is that the matrix elements of 
the effective vertex v^p are in general not necessarily anti- 
symmetric for these energy functionals. Also, for Skyrme 
functionals, one almost always chooses different vertices 
in the particle-hole {vPp^,^) and particle-particle (v'^'^^p) 
channels, and exploits broken antisymmetry of vPP^^ to 
obtain a more versatile effective interaction, for exam- 
ple in the spin-orbit [i^, or spin-spin parts [i^ . The 
situation is similar for the functionals by Fayans et al. 

By contrast, the philosophy of the Gogny force is to 
use the same antisymmetrizcd density-dependent vertex 
anywhere, although in actual calculations terms that are 
very small in SR calculations and at the same time very 
time-consuming to evaluate are often omitted (47j . As 
all standard parameterizations of the Skyrme and Gogny 
interactions use density-dependent vertices, they cannot 
be mapped on a functional that is the strict expectation 
value of a many-body Hamiltonian Almost all rela- 
tivistic mean-field models that are widely used in the lit- 
erature are explicitly set up as Hartree approaches P, ll^l 
without any explicit exchange terms at all, using phe- 
nomenologieal density dependencies and non-relativistic 
pairing energy functionals. 

Note that any local or non-local energy functional that 
contains only terms proportional to integer powers of the 
density matrices can be put into the form of Eq. ([9]) plus 
similar higher-order terms. For the rest of this paper, 
however, we will assume idealized energy functionals that 
are linear and bilinear in the density matrix of a given 
isospin projection, and possibly trilinear with the two 
isospin projections necessarily involved. We postpone the 
discussion of functionals with non-integer powers of the 
density matrices to Paper HI. 

We will not assume antisymmetry of vPp in the for- 
mal manipulations throughout the paper. Owing to the 
intrinsic antisymmetry of k, however, only the antisym- 
metric part of the vertex is probed in the last term of 
Eq. ([9]) and one can always take v'^'^ to be antisymmet- 
ric, which we do here. The results based on a strict HFB 
method can always be easily recovered from those derived 
for a more general bilinear functional simply by enforcing 
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the antisymmetry of v^p and by taking v'''' 



III. PARTICLE NUMBER RESTORATION 

In order to restore good particle number and in- 
clude the correlations associated with the corresponding 
Nambu-Goldstonc mode, it is necessary to extend the 
EDF framework to a multi-reference formalism. This ex- 
tension requires the explicit treatment of the fluctuations 
of the gauge angle of the gap field. This is particularly 
crucial for situations where the symmetry breaking is 
weak or even absent at the SR level, as it is the case for in- 
stance around closed shells or at high spin. The Variation 
After Projection (VAP) method [ll, ^ IH, HO, [HI, [H 
is superior in that respect to the Projection After Vari- 
ation (PAV) one since the latter cannot compensate for 
the spurious sharp phase transition occurrin g at the SR 
level in the weak symmetry-breaking regime jisl. IsgI. [50|. 
[Sll . [s^ l . An intermediate treatment consists of perform- 
ing a projection after a SR-|-Lipkin-Nogami (HFBLN) 
calculation [13, [H^l ■ This corrects for the princi- 

pal defect of the PAV method as it guarantees the pres- 
ence of pairing correlations in the SR state in the weak- 
pairing regime. However, some doubts have been raised 
in the literature about the quantitative reliability of this 
method [13, [111 • The MR calculation could be extended 
further to incorporate dynamical pairing correlations as- 
sociated with fluctuations of the magnitude of an order 
parameter that quantifies the amount of pairing correla- 
tions present in the SR state [H, [H, [H, [5^] . 

An operator that projects out an eigenstate of the par- 
ticle number operator N with an eigenvalue from any 
many-body wave function is provided by [s^ 



1 



2ir 



(10) 



For the purpose of the present paper, it is sufficient to 
consider the simple case of particle number restoration af- 
ter variation. For the sake of transparent notation we dis- 
cuss the formal framework assuming one type of particles 
only. The extension to two types of particles is straight- 
forward and will be mentioned only whenever necessary. 
A normalized projected HFB state is given by 



1*^) = 



(11) 



denotes the overlap of a gauge-space rotated state with 
the unrotated one. The integration interval in Eq. 
can be reduced to [0, tt] using symmetries of the integral 
whenever the SR state \^ui) has a good number parity 
quantum number [s^, [39|, [Hg] . 



A. Energy in the strict PNP-HFB approach 

In the strict PNP-HFB method, the energy is calcu- 
lated as the expectation value of the Hamilton operator 
in the normalized projected state Y^^) 



27r p~iipN 

27rcjv 



(14) 

where we have used that H and N commute and that 
is a projector = p^ _ The energy kernel E[lp\ 

can be easily evaluated with the help of the Generalized 
Wick Theorem (GWT) [13, [sil, which in the canonical 
basis of |$o} gives 



2 
1 



(15) 



In this expression, the normal and anomalous transition 
density matrices between the ket 1$;^) and the bra ($o| 
are defined as 















('&o|<f>^) 




($o|aJ.4l*v> 





(16) 



(17) 



5,-^ ■ (18) 



where the real and positive cjy = ("I'ol^' ) that reads 



The functional kernel E[(p] defined by Eq. ([T5]) has the 
exact same form as the strict HFB energy functional 
E[p, K, K*] given by Eq. ^ except that the SR density 
matrix and pairing tensor (l5][7]) have been replaced by 
the transition ones (|16m8p . Also, the HFB functional is 
recovered from Eq. psp for (p = 0, which amounts to con- 
necting the SR energy and MR energy kernels through 
E[0] E[p,K,K*]. 



1 



($o|i^''|<fo) = ^ d^e-»^^($o|$^)(12) 



provides the weight of the normalized projected state in 
the normalized SR state it is projected from, whereas 



(^oi^-^) = n ' 

Ai>0 



) 



(13) 



B. Energy in the PNR energy functional approach 

Difficulties arise when trying to construct the multi- 
reference energy kernel £[ip] within a true functional 
framework and connect it to the single-reference one. At 
present, there is no ab-initio formalism to derive MR en- 
ergy functional kernels, of which the SR functional would 
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be a special case, and one can only reverse engineer the 
procedure and extend the SR energy density functional 
to the MR level by analogy with the strict Hamiltonian 
case. Based on the strict HFB and PNP-HFB methods 
described above, EDF practitioners have used a proce- 
dure where SYp] = f [p""^, k"*^, *] is postulated to be 
the MR e nergy k ernel that corresponds to a given SR 
functional [ll lii, HO, IMI, HI . In this case, the MR en- 
ergy corresponding to particle number restoration takes 
the form 



?N 



dip 



-iipN 



(19) 



where £[ip] denotes the set of MR energy functional ker- 
nels associated with each gauge angle tp. A kernel £[(p] is 
a functional of the bra ($o| and of the ket \^^), in such 
a way that £^ depends only implicitly on the projected 
state [sgj and cannot be factorized into a form similar to 
the left-hand side of Eq. We will call this procedure 
the "use of the GWT" below, although strictly speaking 
it is not the GWT that is applied, but a formal analogy 
to the extension at play in the strict Hamiltonian case 
when using the GWT that is exploited. 

On the one hand, the standard strategy based on the 
GWT analogy to define the non-diagonal functional en- 
ergy kernel £[ip] from the single-reference functional re- 
placing SR density matrices by the transition ones guar- 
antees that the MR energy functional passes all consis- 
tency requirements thought of so far |2l| . On the other 
hand, this procedure is also at the origin of the diver- 
gences and finite steps discussed in Ref. [H, [1^ . In Pa- 
per I we proposed the general formalism appropriate for 
a remedy of these problems. The remedy is valid for any 
type of multi-reference calculation but is limited to EDFs 
depending on integer powers of the density matrices as is 
further elaborated on in Paper III. The goal of the follow- 
ing sections is to discuss the origin of the problem further 
and to illustrate the general regularization procedure in 
its application to PNR. 

We note in passing that in PNR and all other MR- 
EDF calculations the energy is the only observable that 
is currently determined from a functional; all other ob- 
servables that are routinely calculated within such an 
approach are obtained as matrix elements of the corre- 
sponding operator between projected states, such that 
they do not contain spurious contributions. 



IV. SELF-INTERACTION AND SELF-PAIRING 
A. Single-Reference level 

1. Self-interaction 

Microscopic methods for low-energy nuclear structure 
physics usually describe a self-bound nucleus in terms of 
nucleons characterized by their experimental mass. In 



such an approach, a nucleon should not gain energy by 
interacting with itself. Its so-called self-interaction en- 
ergy, which can be extracted from the one-orbital limit 
of the interaction part of the energy functional £f^ = 
f [p^^,0,0] in the canonical basis, should be strictly 
zero. This requirement is, however, not fulfilled for most 
functionals used in electronic DFT d, i, [H, HO, ^ ^ 
or nuclear EDF methods [lOl- Energy functionals with 
higher-order density dependencies than those discussed 
here might also exhibit multi-particle self-energies, not 
having the proper n-particle limit of the energy func- 
tional [g^ - 

Let us consider the energy of a single Fermion occu- 
pying the canonical state 0^, divided by the probability 
Pflji ~ ^^^^ state to be occupied in the auxiliary 

state |$o) 



2 M/^MM A* ' 



(20) 



This expression shows that a self-interaction arises when- 
ever the vertex v^p is not antisymmetric, vP^i^^ 0, 
which is impossible when calculating the exact matrix 
element of a Hamilton operator, but happens for gen- 
eral energy density functionals. The total one-body self- 
interaction energy is obtained summing all individual 
contributions £„. 



2. Self-pairing 

Beyond the well-known problem of spurious self- 
interactions, there exists a similar problem of spurious 
self-pairing processes which may arise whenever superflu- 
idity is incorporated into an energy functional in a DFT 
or EDF framework. The rationale behind it is that two 
Fermions occupying a pair of conjugated states should 
not gain extra binding through the pairing interaction by 
scattering onto themselves. This requirement constrains 
the two-particle limit of the theory and the contribution 
of a conjugated pair to the many-body energy. To the 
best of our knowledge, the possibility of self-pairing has 
never been addressed before. 

Self-pairing can be easily identified when isolating the 
energy of two Fermions occupying a pair of conjugated 
states {(t>in4'ji} in the canonical basis. We define the 
direct interaction energy of such a pair by removing 
the one-body contributions defined through Eq. ([^0]) to 
£tifi = ^ [{p^^j Ppp } J {'^J?^! '^^^I I {^/?^ I '^p^ }] and 
by dividing the result by the probability P*j to occupy 
the pair in the auxiliary state j^o) 



f f - f- 



- (v''''- - -^v-^ - ) V 



-pp 



The probability P*^ to occupy the pair 



(21) 



(22) 
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is equal to the probabihty of each state to be occu- 
pied, which is a particularity of fully paired quasipar- 
ticle vacua, Eq. H]). In the strict HFB case where 
= ^MApp = "^Wmm' the two terms on the 

r.h.s. of Eq. ((2T|l combine into 



(23) 



using M?, + vf, = 1. The same result is obtained in a 
strict HF method without explicit treatment of pairing 
correlations. The equality of the two-body interaction 
energy (|23p in the HF and HFB case means that a con- 
jugated pair of states {^, p,} does not gain extra direct 
binding by scattering onto itself. Genuine pairing cor- 
relations originate from scattering to different pairs of 
conjugated states and back. 

For most of the standard SR energy density function- 
als used for nuclear structure calculations, however, the 
three terms in Eq. (|2ip can in general not be recom- 
bined into a single one because the vertices entering £pp 
and iS'*" are not related, either by construction or due to 
approximations. Consequently, the direct interaction en- 
ergy of the conjugated pair is not equal to its zero-pairing 
limit as it should be, which gives rise to a spurious self- 
pairing interaction where one has a contribution to the 
energy functional from the scattering of a pair of conju- 
gated states onto itself. 



3. Further discussion 

In a composite system consisting of two particle species 
such as atomic nuclei, the like-particle self-interaction for 
a given particle species is obtained as the one-particle 
limit of the interaction energy for this particle species, 
while keeping the particle number of the other parti- 
cle species unchanged. Otherwise self-interactions in the 
terms that couple the two particle species will be missed. 

The existence of spurious self-interactions was first 
recognized in Kohn-Sham DFT for electronic systems 
[2^. In this context, the construction of self-interaction- 
free functionals has been studied in some detail, see 
Refs. H, [2^, [gO, [m, [ill and references given therein. It 
turns out to be not trivial at all knowing that the stan- 
dard correction method is formulated within the frame 
of so-called orbital-dependent energy density function- 
als [63I . [g^ I and significantly complexifies the calcula- 
tions through the modification of both the total energy 
and the single-particle equations of motion. The (un- 
known) exact Hohenberg-Kohn functional of DFT is of 
course self-interaction free. The spurious terms arise 
when constructing approximate energy functionals that 
are tractable for the use in actual calculations; i.e. self- 
interaction is one of the prices to pay for replacing the 
exact many-body problem by a much simpler set of cou- 
pled one-body problems. It is of course desirable to work 
within a theory that conserves the Pauli principle, but 



its restoration is mandatory only when its violation af- 
fects observables of interest on a scale comparable with 
or larger than the precision desired and reachable within 
a given method. The situation is thus similar to the 
necessity to restore other broken symmetries. As a mat- 
ter of fact, the merits of self-interaction corrected energy 
functionals for electronic DFT are still debated from a 
phenomenological point of view, as they improve some 
observables, but at the same time degrade others when 
compared to uncorrected functionals; see Ref. [g^l and 
references given therein. 

The same remarks apply to self-pairing. Both self- 
interaction and self-pairing processes are actually rooted 
in a violation of the Pauli principle at the level of the 
two-body (or even higher-order) density matrix in the 
definition of the energy functional. It is important to 
stress that they are solely a shortcoming of common en- 
ergy functionals and not of the auxiliary states of refer- 
ence used, as the latter are set up as antisymmetrized 
product states. In particular, all observables other than 
the energy, which are customarily calculated as expecta- 
tion values of the corresponding operators, do not exhibit 
any explicit spurious contributions, although they might 
be indirectly affected through the use of density matri- 
ces that are determined from the solution of a variational 
equation that uses an energy functional containing spu- 
rious contributions as an input. 

In the nuclear context, the possible contamination 
of nuclear energy density functionals by spurious self- 
energies has been noticed before [l], [s^, H^, El] , but was 
never studied in quantitative detail so far. 

It has to be stressed that using self-interaction and self- 
pairing free energy functionals is not per se equivalent to 
the use of an effective Hamilton operator. Indeed, self- 
interaction, as usually characterized, and self-pairing, as 
presently defined, probe only the exchange symmetry of 
a particle in the canonical basis with itself and its con- 
jugate partner, not the exchange symmetry between all 
particles. Asking for a full restoration of the Pauli prin- 
ciple necessarily leads to using a genuine Hamilton oper- 
ator 13011. 



B. Multi-Reference level 

The appearance of self- interaction and self-pairing pro- 
cesses persists to MR calculations whereas new spurious 
contributions particular to the MR level arise from the 
construction of non-diagonal energy kernels. The ex- 
tension of the self-interaction and self-pairing concepts 
to the multi-reference framework, however, is not at all 
straightforward. For instance the very notion of "occu- 
pied" orbitals is ill-defined for transition density matri- 
ces between arbitrary quasiparticle vacua. In the case of 
particle-number restoration, the situation is significantly 
simplified owing to the fact that all vacua entering the 
PNR energy (fTO|) share the same canonical single-particle 
basis, which consequently also is the canonical basis of 
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the Bogoliubov transformation linking any pair of these 
vacua. As demonstrated in Paper I it is precisely the 
latter canonical basis of the transformation connecting a 
given pair of mixed vacua that must be used to mean- 
ingfully identify self-interaction and self-pairing contribu- 
tions to the corresponding multi-reference energy kernel. 



bital 4>^ divided by the probability p*^ to occupy that 
orbital in the projected state \^^) is given by 



1. "Naive" extension of self-interaction 



In the context of PNR multi-rcfcrcncc calculations, the 
energy of a single Fermion occupying the canonical or- 



£N 1.2-n -i<pN ^4 ^iiy:, 

^ = W + k — / 2X22.^ n + ^''^) • (24) 



T ^ 

The one-body density matrix p of the projected state 



(W" V"> C'^^ p-^V" C^^ p-^V" 



is diagonal in the canonical basis of the HFB state it is 
projected from, which means that the canonical basis of 
the underlying HFB state is also the natural basis of the 
projected one. 

As for the SR case, the energy reduces to kinetic 
energy when antisymmetric vertices w'''' are used. How- 
ever, an important aspect specific to the MR case is that 
the integrand appearing in Eq. (|24|) contains a potential 
(simple) pole for Lp = 7r/2 and vf^ = uf^ = 1/2, i.e. when 
the state p, is located at the Fermi level and is not more 
than twofold degenerate in terms of occupation numbers 



v^. If the states present a higher degree of degeneracy, 
an additional factor in the norm overlap will cancel out 
the dangerous denominator. 



S. "Naive" extension of self-pairing 

In multi-reference EDF calculations, the direct inter- 
action energy of a conjugated pair as defined above takes 
the form 



MM 



N 



27r 



dip 



-iipN 



M 



P-P-f^f^ M 



y2 ^2iip 



vle"^^), (26) 



where 



P _ = 

MM 



{■^^\alala-,a,\^^) 



^MM 



(27) 



is the occupation probability of the pair (/i, /i) in the pro- 
jected HFB state. The probability P^^ is equal to the 

probability p*^ of each state to be occupied as we assume 
the underlying SR state to be a fully-paired quasiparticlc 
vacuum with even number parity. 



Using a genuine Hamilton operator, for which w^^^^ = 
^MMMM = ^mPmm = ""t^P^i-^t^ matrix elements entering 
Eq. (|26|) can be recombined in such a way that the poten- 
tial pole disappears [Tsl and that the zero-pairing limit 
is again recovered 



^MM~-^M ~^M 
MM 



TV 



(28) 
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In the EDF formalism, however, the recombination of 
terms in Eq. (pS]) that gives Eq. (|28p cannot be achieved 
anymore. In this case, the integrand in Eq. (|26p contains 
the same kind of pole as the integrand in Eq. ([M]) . 

C. Poles versus "true" self-interaction and 
self-pairing 

In the previous section, we have shown how the self- 
interaction and self-pairing persist to the multi-reference 
EDF framework in the case of particle-number restora- 
tion. What cannot be deduced from such an exten- 
sion of the single-reference case, Eqns. ^U\i and (PT|) . to 
the multi- reference case, Eqns. p4p and (|26p. is if self- 
interaction and self-pairing processes are actually respon- 
sible for the poles. Indeed, recalling our general analysis 
of possible spurious terms in MR energy density func- 
tionals from Paper I, there are in fact two distinct levels 
of spuriosity contained in Eqns. ((M]) and which arc 
of different origins. 



The first level is a consequence of using effective ver- 
tices that are not antisymmetrized, and/or that are dif- 
ferent on the particle-hole and particle-particle channels. 
In the MR framework, such spurious contributions ap- 
pear in the diagonal energy kernels, which are equivalent 
to the self-interaction and self-pairing contributions to 
the SR energy density functional discussed in Sec. IIV Al 
and also enter the off-diagonal kernels. Neither contain 
poles; hence they cannot be at the origin of the diver- 
gences and steps which are the target of the present work. 

In addition to that, a second level of spuriousity arises 
as a consequence of constructing non-diagonal energy 
kernels in analogy with the generalized Wick theorem, 
although strictly speaking the GWT applies only to ma- 
trix elements of operators. As a matter of fact, and as 
demonstrated in Paper I, using a SWT-motivated proce- 
dure rather than a GWT-motivated one does not lead to 
the second level of spuriosity. Taking the example of a 
bihnear EDF, the use of the GWT instead of the SWT 
gives an additional contribution of the form 



cN 



CG 



27r ^—iipN 



dip 



(29) 



E\l(vPP +vPP-- +vPP +vPP- )-v'^'^ ] (n V ) 

M>0 



dip 



that is absent in a SWT-motivated procedure and which contains a pole clearly similar to those discussed in connection 
with Eqns. (|24ll26p . Having identified the contribution ((29)) caused by the use of the GWT, we defined in Paper I the 
regularized MR energy and energy kernels, respectively, as 



-TV 



^REG — '-CG ' 

£regW\ = ^cg\hA ■ 

Removing £qq from Eqns. ([M)) and one obtains the "true" MR self-interaction 



(30) 
(31) 



cN _ 
^S/ — 



It: ^-iipN 



dip 



Z7l Ca 



/i>0 



-N 



PP 



2ir ^-iipN 



dip 



27r c% 



vl + v,\^^^) +2uXie'^^ - 1)1 n i^l + -I ^'-n ■ (32) 



i,>a 



and the "true" self-pairing contribution 



cN _ 
^SP — 







2ir ^-iipN 



2ttc% 



[ MMW 2 



M>0 



v>0 



27r ^-iipN 



dp> 



2ttc% 



,2 „,2 



2 / 2 I 2 2iip\ I / 4 2 2 4 \ / 2iip -, \ 



(33) 



both of which belong to the first level of spuriosity and do not contain any dangerous poles. The expressions 
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FIG. 1: (Color online) Particle-number restored deformation 
energy surface of calculated with SLy4 and a density- 
dependent pairing interaction and the corresponding single- 
particle spectra of protons and neutrons as a function of the 
axial quadrupole deformation for i = 5 and 199 discretiza- 
tion points of the integral over the gauge angle (lowest panel) . 
There are clear anomalies that appear when either a proton or 
neutron single-particle level crosses the Fermi energy. The di- 
mensionless quadrupole deformation /32 is defined in Eq. (|66|l . 



([5^ and ([55)1 could also have been obtained directly from 
Eqns. (79) and (80) of Paper I. 



D. Impact of the poles on PNR energies 

In the previous section, we demonstrated that the spu- 
rious contribution contains poles. Figure [T] illus- 
trates, through a realistic calculation of the particle- 
number restored deformation energy surface of ^^O, the 
impact of such poles for a functional containing a frac- 
tional power of the density matrix. The SLy4 parame- 
terization of the standard Skyrmc EDF is used in con- 
nection with a density-dependent pairing energy func- 
tional, which was used in many MR calculations be- 
fore [63, m, m, H [Zll, [zS, [Zi . In practice, the integral 
over the gauge angle appearing in Eq. (jl9p is discrctizcd 
into a sum using the Fomcnko expansion, as will be ex- 
plained in Sec. IVI Bl below. It is important to stress that 



all observables calculated as operator matrix elements, 
e.g. particle number, quadrupole moment, radius, etc., 
are converged using five integration points. The particle- 
number restored energy functional, however, does not 
converge. Instead, one observes the development of sev- 
eral localized divergences as one increases the precision of 
the calculation, which appear exactly where neutron or 
proton levels cross the Fermi energy; i.e. where their oc- 
cupation probability is = 0.5. In spite of the evidence 
for their appearance presented in Refs. [H, [H, [H, [2^, 
the divergences remained undetected so far in our PAV 
calculations, because on the one hand the appearance of 
the divergence requires a number of integration points 
far above the one used in practical calculations, and be- 
yond what is tractable in connection with other projec- 
tions and mixing of different deformations, and because 
on the other hand the divergences are sufficiently local- 
ized in deformation space that the area obviously affected 
by the pathology is smaller than the typical distance of 
states commonly used when calculating energy surfaces 
and when mixing states with different deformations. 

At this point, three questions arise (1) do the diver- 
gences seen in Fig. [T] constitute the only pathological 
manifestation of the poles? (2) Do divergences mani- 
fest for any type of functional, i.e. irrespective of the fact 
that it is bilinear, trilinear or contain non-integer powers 
of the density matrices? (3) Is the spurious contribution 
isolated in Eq. l(29)) responsible for all problems associ- 
ated with the poles; i.e. would removing it from PNR 
energy kernels properly regularize the MR EDF calcu- 
lation? Answering theses questions will be the aim of 
Sec. I VII Before discussing the results obtained using the 
method proposed in Paper I to regularize MR energy 
kernels, we discuss the pathological manifestations of the 
poles in more detail through a complex plane analysis, 
following Ref. 0. 



V. COMPLEX PLANE ANALYSIS 

The integral over the real gauge angle can be reformu- 
lated as a contour integral in the complex plane, which 
allows the analysis of the energy functional in terms of 
its poles within the integration contour [2^. In fact, 
particle-number projection was first introduced through 
such complex contour integrals [13, [zl] ■ It was only af- 
ter Fomenko [tH] demonstrated that a simple trapezoidal 
rule gives a very efficient discretization of integrals over 
the gauge angle that Eq. pT|) became the standard way 
to formulate and evaluate PNR observables. 



A. Analytic continuation 

To that aim, one introduces the complex variable z = 
e"^. As a result, quantities used in the PNR method 
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Rc(z) 



FIG. 2: Schematic view of the analytical structure of the tran- 
sition densities defined in Eqs. (|38II40P and of the PNR func- 
tional energy kernel £[ip] in the complex plane. Poles marked 
with filled circles are within the standard circular integration 
contour of radius R — 1, while those outside are marked with 
open circles. The cross marks the location of the SR energy 
functional at z = 1. The operator e"^^ produces a rotation 
in gauge space, while e''^ is a shift transformation as defined 
in Eq. (|45ll . 



1)^ 
(34) 

(35) 

(36) 



involve an integration over the unit circle Ci (|z| 







dz 1 


/ 

JCi 


2iTTCN z'^+^ 


£^ = 




dz £ [z] 


/ 

JCi 


2iTrc% z^+i 


c% = 


/ 

JCi 


dz 1 
2^7^ z^+i ^ 



whereas the overlap now reads 



M>0 



(37) 



Finally, the transition density matrix and pairing tensor 
extended to the complex plane become 











Z2 








ul + vl 




^,zl * 

'^fj.f — 


u^v^ z 


2 




Z^ 



5uL 



(38) 
(39) 
(40) 



B. Energy functional kernels 

Taking advantage of the Cauchy residue theorem, go- 
ing to the complex plane allows the calculation of all 
quantities of interest in terms of poles of the integrand 
located inside the integration contour. For the norm 



-TV 



7^es(0) 



M>0 



(41) 



or any other operator matrix elements between projected 
states, only the pole at z = contributes. 

The situation is different for the PNR energy as addi- 
tional poles at finite z^ = ■izi\u^\/\v^\ enter the energy 
kernel £ [z]. Thus, Eq. ([55]) takes the form 



E 



— TZes{zi) 



Zi=0,|z±|<l 



'N 



£[z] 

■yN- 



(42) 



with contributions from the pole at the origin and from 
all pairs of "hole- like" poles at z^. The situation is 
schematically depicted on Fig. [2] The location of the 
pole associated to a given pair (/i, fi) moves along the 
imaginary axis as the occupation changes with defor- 
mation. When the corresponding pole crosses the unit 
circle, either entering or leaving the Fermi sea, the in- 
tegrand is non-analytical on the integration contour and 
the integral diverges. 

The point has now come to realize that the divergences 
constitute the most obvious part of the problem, but do 
not contain the entire problem. As can be seen from 
Eq. (|42)) . the poles at |z^| < 1 contribute to the en- 
ergy when using an energy functional that contains self- 
interactions and self-pairing. On the other hand, only 
the pole at the origin contributes in the strict PNP- 
HFB/Hamiltonian framework as the poles at \z^ \ do not 
exist in this case. Consequently, one has to ask the ques- 
tion whether or not the contributions from the poles at 
< 1 2^ I < 1 to the projected energy are physical, in par- 
ticular when realizing that the contribution of a given 
pole can be many orders of magnitude larger than the 
total energy gain from PNR [2^. In addition, a pole 
at finite \z^\ entering or leaving the integration circle 
does not only provoke a divergence but also provides the 
PNR eriergy with a finite step after the crossing is com- 
pleted [2^. Looking carefully at the potential energy 
surface obtained using L = 199 integration points, such 
a step can be seen in Fig. [TJ i.e. compare the energy be- 
fore and after the crossings at (32 = +0.22 and (32 = —0.3. 
As a matter of fact, the binding energy jumps from one 
potential energy surface to another. 



^ We abusively replace the gauge angle by the complex variable 
2 in all our expressions; i.e. SR states characterized by the gauge 
angle i/j, are extended into |<I>z) to denote SR states any- 

where on the complex plane. In particular, the unrotated SR 
state, denoted as |3>o) when using </? as a variable, is written as 
It&l) when using z as a more general variable. 



C. Spurious contributions 



In Section IIV C[ we have identified £^q as the only 
possible source of spurious poles. In order to obtain a 
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deeper insight to its content, we rewrite Eq. ((29|) as 

cN 
^CG 

I dz £cG H TT/ 2 , 2 2n 



L2 V^A" 

M>0 



X 



dz z^-1 



(43) 

and define in passing the spurious contribution £cg [z] 
to the MR energy kernel over the entire complex plane. 
From Eq. (|43|) , the spurious contribution of each pole to 
the PNR energy can be calculated. As for the total en- 
ergy, the poles of the integrand are located at zq = and 
z,t = ±i\u,,\l\' 



that removing from £^ 



This has the important consequence 
^fjQ iiuiii does not only extract the 
contribution of the poles at \z^\ < 1 but also a spurious 
contribution of each conjugated pair (/x, /2) to the physi- 
cal pole at Zq = 0. The latter could not have simply been 
guessed from the analysis of the analytical structure of 
£[z\ in the complex plane. As a matter of fact, the spu- 
rious contribution from the pole at zq = is absolutely 
essential for the internal consistency of £cg- '-"^'^ 
hand, it was shown in Ref. [1^ that the energy associ- 
ated with a single pole at |zj^| < 1 can be gigantic (away 
from where it might be divergent). On the other hand, 
the total spurious energy hidden in a PNR method can- 
not be larger than the energy gain from particle number 
restoration itself, which is on the order of at most a few 
MeV. It is only the combined contribution from the poles 
at Zq = and z^, which nearly cancel each other, that 
will give reasonable values to the total spurious energy 
£qq as will be exemplified below. 



The residue for the pair of poles at 
Eq. ([33]) can be evaluated analytically 

ne^G{4) 



TZes{zi) 



■>0 {ul 



contained in 



,,2 ^2N 



Af+2 



l + (-l) 



N 



2i 



N 



n 



2 2 



9 9 



(44) 



Note that TZe^Q{z^) is zero if projecting on an odd par- 
ticle number TV as the underlying reference state ([T]) has 
been chosen to have an even number-parity quantum 
number [H, [sll . The generalization of the present dis- 
cussion to the case one- (or 2n -|- 1) quasiparticle states 
with an odd number-parity is straightforward, but not 
important for the purpose of this paper. 

The total contribution from the pair of poles < 
|z^| < 1 to the PNR energy is then obtained by re- 



placing the integral in Eq. P5|) by 2i7r TZe^Q 



^), where 



TZc^qIz^) is given by Eq. (|44)) . We will discuss the in- 
dividual contributions from the poles in Sec. IVII below. 
Note that calculating the residue of the pole at zq is much 
more involved because it is a pole of order -t- 1. Its 
residue can in fact be calculated analytically through a 
recursive formula, which, however, involves a sum over 
such a large number of terms that it is of no practical 
use and is not reported here. In any case, one can access 
the spurious contribution from the pole zq by subtract- 
ing the analytic expression of Eq. (|44p from a numerical 
evaluation of the full expression given by Eq. (j^. 



D. Properties under shift transformation 

The interpretation of the poles at z^ 7^ becomes 
clearer when looking at the properties of the PNR energy 
functional under a so-called shift transformation p5| . In 
the present paper, we choose a slightly different definition 
of the shift transformation from the one used in Ref. 123 



(45) 



such that the shift transformation operator 6^''+*'''^^ used 
in [25I is the product of ours (j45p and a rotation in gauge 
spacc^. In contrast to a gauge-space rotation that is 



unitary, the shift transformation (j45p is non-unitary and 
changes the norm of the product state. 

In the complex plane, the shift transformation (j45p 
corresponds to a radial shift of z from z ~ e*'^ to z' = 

e'^'f , see Fig. [51 Thus, projecting a shifted HFB state 
on particle number amounts to changing the radius of 
the integration circle from R ~ 1 to R ^ j25| 



ip-iril 



1 



dz 



R 



N 



2m z^+ 



(46) 



where we have made the substitution z' — e^'^ in the first 
line and the substitution z = Rz' in the second one. Both 
expressions will turn out to be useful below. The overlap 
between the non-normalized projected SR state and its 
counterpart shifted along the real axis is given by 



{R)^{^,\P''\<^R 



(47) 



with c\i as defined through Eq. ([T2|) : i.e. c 



(!)■ 



All normalized projected matrix elements are shift in- 
variant if the operator O in question commutes with N. 
Just as the exact ground-state energy, its approximation 



^ Starting from a circular contour, the additional rotation in the 
definition of Ref. [Ssl does not make any difference. The situation 
would have been different if we had started from a non-circular 
contour. 
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obtained through the particle-number restored expecta- 
tion value of the Hamilton operator is shift invariant. On 
the other hand, this is not the case for standard particle- 
number restored energy density functionals (25j . The vio- 
lation of shift invariance is obviously a consequence of the 
presence of the poles at finite contained in the PNR 
energy kernel constructed on the basis of the GWT. For 
a given spectrum of poles the energy £^ changes by 
a finite quantity whenever the integration circle crosses a 
pair of poles \z^ \ in the course of a shift transformation. 
As a result, the PNR-EDF is shift invariant only over a 
finite range of values of the shift parameter r] [25[ . This 
result clearly points to the unphysical nature of these 
poles. 



E. Sum rules 

One might wonder where the energy that is 
added/removed when crossing a pole with the integra- 
tion contour comes from/goes to. In the present section, 
two different sum rules involving PNR energies £^ ex- 
tracted from a given SR functional are carefully derived 
and discussed to answer such a question. 



1. Radius-weighted sum rule 

As it is introduced in Ref. [1^, we first discuss 
the characteristics of the radius-weighted sum rule 
^ cj^iR) £^ {R), although we already insist here that the 
physical sum rule of interest is the non-radius-weighted 
one discussed in Sec. IVE 21 below. The number R ap- 
pearing in the sum rule is taken to be real even though it 
is possible to formulate the sum rule using an arbitrary 
complex number of norm R [2^ . Our conclusions will be 
insensitive to this detail. 

First, let us recall how such sum rules arise in the 
operator- and wave-function-based context. Inserting 
the complete set of normalized particle-number projected 
states'^ 

^|*^)(.|/^|=^P^ = 1 (48) 

N>0 N>0 

into an unprojected shifted matrix element of an operator 



^ The fact that one does not need to sum over N < can be seen 
as a consequence of the fact that [^^) = for A'^ < as a result 
of the disappearance of the physical pole at 2 = in the contour 
integral of Eq. I|34|l . Note that the normalized projected-state on 
AT = is \^°) = |0>. 



O that commutes with N gives 

{^i\6\^r) = (<i>i|6e''*|<I'i) 

= ^($i|Oe''*|*^)(*^|$i) 

JV>0 

= ^c^(i?)0^, (49) 

N>0 

where we have used that e^^l^i/^) = R^\'i>^) and define 
= Equation ^ expands the 

shifted SR matrix element 0[R] = ($i|0|$fl) in terms of 
average values of the operator in all normalized pro- 
jected states. Applied to the Hamilton operator, Eq. ([49]) 
reads 

E[R]=J2cUR)E'', (50) 

Af>0 

and provides for 77 = (i? = 1) that the strict HFB 
energy decomposes into strict PNP-HFB energies (with 
iV > 0) weighted by the probability to find the normal- 
ized projected states into the SR state. In Eq. (|50p . 
the sum could be further reduced to > as the 
contribution from the term A = is strictly zero, i.e. 
Cg = E[z = 0] niy>o ^ ^- Such a result rehes on 
the fact that only the physical pole at z = contributes 
to the integral providing E'^ . 

Let us now come to the EDF context and lay out some 
specificities that are crucial to provide a meaningful dis- 
cussion of sum rules, (i) In Eq. (|50p . it was not neces- 
sary to specify the integration contour used to calculate 
E'^ as the latter is shift invariant. In the EDF context 
where the shift invariance might be broken, it is manda- 
tory to specify the contour employed. Consequently, the 
notation 8^ (R) is used whenever necessary to character- 
ize that a circular contour Cr of radius R is employed 
to calculate PNR energies, (ii) There is no equivalent 
to "inserting a complete set of states" in the EDF con- 
text as one directly postulates the PNR energy under the 
form of a functional built from one-body transition den- 
sity matrices and integrated over the gauge angle, and 
not from the expectation value of a Hamilton operator in 
projected many-body wave functions. As a consequence, 
the existence of a sum rule similar to the one discussed for 
operators is neither obvious nor trivial. By contrast to 
the above derivation, one has to start from the weighted 
sum over PNR energies and see if and how it recombines 
in the same manner as for an operator matrix element. 
To obey a sum rule analogous to the one provided by 
Eq. (|50[) can thus be demanded as a consistency require- 
ment for MR energy density functionals. To recover the 
SR energy from such a sum rule, it is a necessary condi- 
tion (but not sufficient) that the MR energy kernel £[z] 
is set up such that it gives back the SR energy functional 
£[p,K,K*] for z = 1, as assumed throughout this paper, 
(iii) The sum rule considered in the present section actu- 
ally differs from the one discussed in Ref. [2^. Indeed, it 
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is mandatory in the EDF context to make the sum run- 
ning over both positive and negative "particle numbers" . 
As will be shown below, the latter is crucial to establish 
the expected sum rule when individual particle-number 
restored energies £^ are not shift invariant, i.e. when MR 
energy kernels £[z] possess spurious poles at finite . In- 
deed, the product c%{R)£^{R) is different from zero in 
this case for < because, although the physical pole 
at z = disappears from the integrand as it should, the 
poles at finite contribute. This is certainly the most 
direct proof of the non-physical nature of such poles and 
non-regularized energy functionals. In the context of the 
real-space derivation of Ref . [2^ , obtaining the appropri- 
ate sum rule calls for using the correct Fourier decompo- 
sition of the periodic delta function over all irreducible 
representations of C/(l) including those characterized by 
negative integers N; i.e., J2n=-oo e"**^^ = 2TrS2TT{v)- In 
the following we proceed in the complex plane to estab- 
lish the needed sum rules. 

First, the change of variable z — Rz' is performed in 
order to recover an integration over the unit circle 



+°° J- £[Rz] context as a matter oi tact, tne contrioution irom tne 

^ c%{R)£'^{R)^ J2 f poles of at z± depends on the original contour Cr 

— M— — <^i^'^ flTiH nn the infinitpsim^j.l shift trpnsfnrmpitinns lej^.Hinp" tn 



The physical pole at z = 0, which is of order N + 1 
in £^ , has transformed into two simple poles at z = 
and z = 1 in both integrals in Eq. ((52|) . Note in passing 
that the pole at z = would have not appeared if we had 
grouped the component iV = to the sum over positive 
numbers. The pole at z = 1 is on the unit circle and is 
thus located inside of Ci+, but outside of Ci- . Thus, it 
contributes to the first integral only in Eq. (|52p and pro- 
vides the sum rule with the contribution £[R] {^i\^h) 
which represents the transition kernel involving the orig- 
inal HFB state |$i) and the state \^r) shifted along the 
real axis to z = i?. 

In the strict PNP-HFB method, this is the only con- 
tribution to Eq. ([52]) as the residue of the simple pole 
at z = 0, which corresponds to the contribution from the 
N — component, is zero for the reason explained earlier. 
In any case, such a pole contributes to both integrals in 
Eq. ([52|) such that any finite residue would have canceled 
out anyway. Thus, the sum rule (jSOp is recovered. 

The question is whether this still holds in the EDF 
context As a matter of fact, the contribution from the 



N=-CyO 



N=-ca 



(51) 

We recall that £^{R) is proportional to l/cj^(R), 
Eq. (fT9|) . As a consequence, c^(i?) = alone is not a 
sufficient condition that the contribution of a given A^ 
to the l.h.s. of Eq. ((5T|) vanishes, as c%{R)£'^ (R) might 
remain finite. We will come back to this below. 

To invert the summation and the integral in Eq. (|5ip 
and perform the summation explicitly, the power series 
must be (uniformly) converging on the integration con- 
tour. To ensure this property, one has to separate the 
sums over positive and negative N and use the (local) 
shift invariance of £^ to scale the integration radius ap- 
propriately in each of the two terms thus generated. Us- 
ing two infinitesimal shift transformations characterized 
by 77+ > (77- < 0) for A^ > (A^ < 0), the right-hand- 
side of Eq. (|5ip splits into two geometric series converging 
separately and uniformly on the corresponding integra- 
tion contours Ci+ (Ci-). Performing the summation of 
both geometric series, one obtains 



and on the infinitesimal shift transformations leading to 
Eq. ((52)) . If the shift transformations are such that no 
pole appears in between the two contours Ci- and Ci+, 
all poles with |z^| < R contribute to both integrals and 
cancel out in Eq. whereas all poles with |z^| > R 
do not contribute to cither of them. This proves that, 
except for the ill-defined case of a pair of poles sitting on 
the original integration circle Cr, one can always perform 
two infinitesimal shift transformations to prove that 



-t-00 

c%{R)£''{R)^£[R]{<i>^\<^n) 

N=-oo 



(53) 



+ 00 



r{R)£''iR) 



N= 



c 



1+ 



£[Rz] 



2iTT z{z — 1) 



(52) 



Equation ([55]) thus expresses that the expected sum rule 
is found to be valid, even for contaminated and yet un- 
corrected EDFs, i.e. using energy kernels constructed on 
the basis of the GWT, at the price of including the con- 
tributions from unphysical components (A < 0). 

Applying the same derivation as above to the spurious 
contribution isolated in Eq. (|43p . one obtains 
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+00 

^ c%{R)£^a{R) = £cg[R]{^i\^r) (54) 



(R^ \^ (yPP +v''-^--- +V''''- - +vPP- ) -v^"" 1 TT fu^ + i?' 



u>0 



which is zero for i? = 1 as 2; = 1 is the only point in the 
complex plane where the GWT-related spurious contri- 
butions to the MR energy kernel is zero. 

It is crucial to analyze further the cancellation of the 
contribution of spurious poles in Eqs. (|52ll53p . Indeed, 
such a cancellation relies on the original summation over 
both positive and negative " particle numbers" in the def- 
inition of the sum rule. If one sums over positive particle 
numbers only, all pairs of poles situated inside Cr con- 
tribute to the sum rule. This is puzzling as it is clearly 
unphysical to consider negative " particle numbers" . In- 
deed, one necessarily has c1j{R) (R) = for < 
when employing a genuine Hamiltonian. However, the 
product cj^{R) £^ {R) is different from zero for < 
if £[z] possesses poles at finite j^^j < R. This is to our 
opinion the most direct way of stating the non-physical 
nature of those poles. In any case, and as proven above, 
one can at least recover a sum rule for uncorrected func- 
tionals at the price of summing over both positive and 
negative particle numbers. If summing over positive val- 
ues only, one obtains, using our example of a bilinear 
functional, 



oral orders of magnitude smaller than the contribution 
from iV > in realistic cases, such that it might pass as 
numerical noise to the unsuspecting eye. 

Subtracting Eq. ([54]) from ((53)) provides the quan- 
tity E 



(R) [£^{R) - £SGiR)] by which the sum 
rule is modified when regularizing the MR energy ker- 
nels. One observes that the non-physical components 
are zero, i.e. cj^{R) £^^q{R) = for < 0, and that 
the sum rule matches the regularized kernel at z = i? 
£reg[R] 



2. Non-radius-weighted sum rule 

The sum rule (|53p is of particular interest when the 
unit circle Ci is used as an integration contour to define 
PNR energies. Indeed, Eq. ([55)) reduces in this case to 



+00 

E 

Af=-oo 



2 cN / r> 
C^t (R- 



-1) = £[z = l]= £[p,k,k], (56) 



^4(i?)5^(i?)-5[i?](<f>i|<f>«) 



N>0 



J2 ^es(z±/i?) 



£[R'. 



z{z-l) 



M>0 



r.PP 



, V'-'--- -\- V'- - -\- V'-' - I — V^- - 
p^ppp ' PP'P'P' PPPP P'PP'PJ PPPP-i 



Ul R'^ vl -pr 

•il + R'^vl J-i 
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which shows that the physical sum rule {N > 0) is broken 
by a finite amount that relates directly to the presence 
of spurious poles at finite z^ inside the original integra- 
tion circle Cr. Note again that the simple pole at z = 
docs not contribute as its residue is zero. Equation ()55[) 
proves that the sum rule derived in Rcf. [25| is incor- 
rect for the cases of interest. In particular, computing 
Eq. (|55p for i? = 1 provides the non-zero amount by 
which the decomposition of the SR EDF into its physi- 
cal PNR components {N > 0) is broken, already for the 
standard integration circle. However, as we will show in 
Sec. IVI D 41 below, the contribution from iV < is sev- 



which expresses that the SR EDF decomposes into PNR 
energies obtained for all possible "particle numbers" 
^ 0. This decomposition actually relics on the (re- 
quired) connection between the SR EDF and the MR en- 
ergy functional kernel; i.e. £[z = 1] = £[p, k, k*]. Equa- 
tion (|56)) is valid prior to any regularization of the PNR 
energy kernel, as long as the sum runs over both positive 
and negative particle numbers. The null sum rule ()54p 
at i? = 1 shows that regularizing the PNR EDF method 
(^55) through the removal of £^q from £^ consists, for this ra- 
dius, of rcshuffiing contributions among different particle- 
number restored energies, in such a way that the decom- 
position of the SR EDF into its physical PNR compo- 
nents (A^ > 0) is fulfilled. Note that the regularized sum 
rule matches the SR EDF precisely because the regular- 
ization does not modify the energy kernel £[z] for z = 1. 

Still, the radius- weighted sum rule considered in 
Sec. IVE2I and in Ref . [25| docs not allow us to study the 
shift invariance of Eq. (|56p . which is the real question 
of interest. Indeed, what matters is whether or not the 
standard decomposition of the SR EDF into c^-weighted 
PNR energies is valid independently on the radius of inte- 
gration chosen initially to compute £^ . In a Hamiltonian 
and wave function based framework, such an invariance 
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reflects the trivial identity 



($i|i7|<i>i) = J2 



N\i2 



N>0 



JV>0 



Translated to the functional framework, this amounts to 
considering the non-radius- weighted sum rule 

(58) 

where c^(l) = cj^ and where the circle of integration 
Cr is the one chosen to calculate PNR energies. Again, 
the power series must be split into two parts to perform 
the summation over particle numbers explicitly. The ini- 
tial circle of integration Cr being above/below the unit 



circle, one needs to perform a finite shift transforma- 
tion to bring the circle associated with negative/positive 
particle numbers on the other side of the unit circle, in 
order to make the corresponding series convergent. If 
particle- number restored energies are shift invariant, one 
can proceed without any difficulty and obtain the trivial 



result that the sum rule X]7V= 



-oo '^N 



is valid independently on the original radius R. This is 
of course the case for a Hamiltonian- and wave- function- 
based PNR method which, once again, would only require 
the summation over positive particle numbers in the first 
place. 

Of course, problems arise if particle-number restored 
energies are not invariant as the shifted circle crosses a 
spurious pole at z^, i.e. if there are poles located in 
between Cr and Ci. Indeed, proceeding to the required 
shift transformation brings an extra contribution to the 
sum rule in this case. Exemplifying the problem for a bi- 
linear functional and an initial radius R > 1, one obtains 



J 



+ 00 



£[p,K,K*]+ c%£^a{R), 



N=-oc l<c\zi\<R 



£[z] 



Ai>0 



(59) 



N= — oo 



with 



N=-CX3 



where Tie^Q (z^ ) is given by Eq. ([44]) and where the sums 
run over all pairs of poles located in between the unit cir- 
cle Ci and the integration circle Cr. Note that, in agree- 
ment with Eq. ([54)) . the sum rule (|60|) is zero for i? = 1 as 
no pole resides between Cr and Ci in this case. However, 
it is easy to see from Eq. pi]) that X]Ar<o '^^cai^^) ^ 
diverging geometric series of common ratios | | > 1 for 
R > 1; i.e. the sum rule is broken by a diverging amount 
as soon as poles are located in between the integration 
circle Cr and the unit circle Ci . One can check that the 
situation is similar if i? < 1 and the conclusion identi- 
cal. Regularizing the PNR EDF through the removal of 
£qq{R) amounts to transferring the second term in the 
right-hand side of Eq. (|59p to the left-hand side. Doing 



so restores the physical value {£[p, k, k*]) and the shift 
invariance of the sum rule as the shift invariance of each 
individual PNR energy £^ (R) is actually restored. As 
c\ £^^q{R) — for < 0, the sum rule is in fact re- 
stored and made shift invariant by summing over positive 
particle numbers only 

^4f^sG(i?)-£:[p,«,K*]. (61) 

Last but not least, it is of interest to look at the non- 
regularized sum rule obtained by summing over physical 
components only (A > 0). In this case, the physical sum 
rule calculated for i? > 1 is broken by a finite amount 
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with 

+00 2 2 2 2 

E ^^(^) = E (^;^^.. + -^^^ + + -z,^ " ^;:;^..] n "'^""7'^"^ ^ (63) 



where the sum runs over all pairs of poles located inside 
the integrations circle Cr. This time, however, and as 
already made clear above, the sum rule (|56p is not even 
recovered for i? = 1 as the last term of Eq. ([62|) does not 
go to zero. Regularizing the PNR EDF through the re- 
moval of amounts to transferring the second term in 
the right-hand side of Eq. ([62|) to the left-hand side. Once 
again, doing so restores the physical value, i.e. £[p, k, k*], 
and the shift invariance of the sum rule. 



3. Main conclusions 

The first conclusion is that the decomposition of the 
SR energy f [p, k, k*] into its physical (N > 0) particle- 
number restored components is (i) always fulfilled for 
a Hamiltonian- and wave-function-based method, what- 
ever the chosen integration circle is, while it is (ii) broken 
by an amount that depends on the chosen integration 
contour for an EDF-based PNR method if MR energy 
kernels £[z] contain poles at finite z^, but (iii) recov- 
ered for any value of R after regularizing £^ through the 
removal of £qq . 

The second conclusion is that the decomposition of 
£[p,K,K*] involving unphysical components {N < 0) is 
(i) always fulfilled in a Hamiltonian- and wave- function- 
based PNR method as unphysical components do not 
contribute anyway (ii) fulfilled in the EDF context if in- 
tegrating over the unit circle Ci, even for MR energy 
kernels £[z] plagued by poles at finite z^ (iii) fulfilled for 
any integration circle Cr by the regularized EDF-based 
PNR method, noticing in addition that unphysical com- 
ponents do not contribute anymore. 



VI. APPLICATIONS 

A. General remarks 

As seen in Sec. |IV] there are two distinct classes of 
spurious contributions to a multi-reference energy den- 
sity functional. The first one represents the "true" self- 
interaction and self-pairing processes which already ap- 
pear at the singlc-rcfcrcncc level. It does not provide MR 



energy kernels with poles; hence, it does not cause diver- 
gences or steps in the PNR energy and does not break 
its shift invariance. The second one is due to the use of 
the GWT out of its context to define MR energy func- 
tional kernels from an underlying SR EDF that contains 
self-interaction and self-pairing contributions. 

As outlined in Sec. IIV A H correcting consistently for 
the standard (true) self-interaction £gj, Eq. ((32|) . is not 
an easy task; the correction enters the variational equa- 
tions already on the single-reference level and leads to a 
state-dependent single-particle field [1^, EO, lU, H^l . The 
same would hold regarding the correction for true spu- 
rious self-pairing £sp, Eq. ([55)) . For that reason, and 
because such spurious contributions are not responsible 
for divergences and steps in the PNR energy, we concen- 
trate here on £qq, Eq. (|29p which is at the origin of the 
specific and dramatic pathologies encountered in PNR 
EDF calculations. Note that subtracting £qq from the 
PNR energy will also modify the variational equations 
of a VAP calculation. Here, we confine ourselves to an 
analysis of the poles and of their impact on the particle- 
number restored energy after the variation. In this case, 
£cQ is easily subtracted a posteriori. 

There is one important limitation to the applicability 
of the regularization method proposed in Paper I and 
applied in the present work. Although it is straight- 
forward to extend Eq. ((29|) to an EDF depending on 
any integer powers of the density matrices, this is not 
the case for EDFs depending on non-integer powers of 
the densities. This is a significant limitation, consider- 
ing that most successful modern functionals use density 
dependencies of non- integer power'*. Indeed, this allows 
them to provide a good description of the most impor- 
tant nuclear matter properties with a very small num- 
ber of terms and coupling constants to be adjusted phe- 
nomcnologically [l|. Also the widely used Slater approx- 
imation to the Coulomb exchange term falls into the cat- 
egory of a density dependent term of non-integer power. 
We analyze the spurious contributions to such category 
of functionals in Paper III, complementing the study of 



* An exception is the relativistic functional [76| used in the MR 
calculations of Niksic et al. |20|| . 
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Dobaczewski et al. (25[. In the present work, however, we 
use instead the particular early parameterization SIII [T^I 
of the Skyrme EDF that contains only bilinear and trilin- 
ear terms in the normal density matrix. We complement 
the SIII energy functional with a density-independent lo- 
cal pairing functional that is bilinear in either the neutron 
or proton anomalous density matrix. For the Coulomb 
energy functional, we only consider the direct term and 
neglect the approximate exchange term that was consid- 
ered in the fit of SIII. As a consequence, all calculated 
nuclei will be underbound by a few MeV, but this is of no 
importance for the purpose of the present paper. Hav- 
ing said that, it is clear that the construction of high- 
precision correctable EDFs, i.e. only containing integer 
powers of the density matrices, represents an important 
task for the future^ . 

The calculation of the various contributions to the cor- 
rection £qq is outlined in Appendix [X] The trilincar 
terms in the SIII functional are motivated by a local zero- 
range three-body force which excludes terms of third or- 
der in the same nucleon density; it only contains terms 
of the kind p'f^{r)pp{r) and Pp{r)pn{r). From a practi- 
cal point of view, the absence of a genuine term of third 
power in the same density matrix has the advantage that 
we do not have to invoke the corresponding correction 
term outlined in Paper I. Instead, the correction of the 
trilinear terms has the structure of the one of bilinear 
terms times the projected density of the other species as 
outlined in Appendix IB 21 

B. Numerical Implementation 

In practice, the integrals over gauge angles arc dis- 
cretized with a simple ??-point trapezoidal formula 

- r d^f{e^^)^lj2f{^'^) (64) 

where we assume the projection of a state with even num- 
ber parity on even particle number to reduce the integra- 
tion interval to [0, tt] . As was shown by Fomenko [t^ , this 
simple scheme eliminates exactly all components from the 
SR state which differ from the desired particle number 
N by up to ±2(L — 1) particles. Although the spread in 
particle number is large compared to the total particle 
number, already small values for L, ranging from 5 in 
light nuclei to 13 in heavy ones, are sufficient to obtain 
a converged projected state. 



^ In practice, one will have to restrict the form to rather low or- 
ders in the density matrices. For example, the EDF recently 
proposed by Baldo et al. l78l includes terms up to fifth power 
in the total density p(r), which clearly lead to self-interaction 
terms |30|| that will require a regularization containing quadru- 
ple sums over single-particle states, which will be too costly in 
realistic calculations. 



It is customary to use an odd number of discretization 
points L in the interval [0, tt] to avoid numerical problems 
that may appear at = 7r/2. This practice docs not 
relate to the real divergences of the energy functional 
contained in £qq that we discuss here, but avoids the 
implicit division of + v'^^ e''^^^ contained in an opera- 
tor kernel by the same factor in the normalization factor 
when evaluating projected operator matrix elements 
(as, for example, particle number, deformation or radii), 
which numerically will not give the analytical result 1 
when u'j^ comes very close to w^. Of course, the numer- 
ical representation of the pole contained in the energy 
functional would not be very precise in this case either. 

With a small modification, the discretization ([64|l can 
also be used to represent complex contour integrals with 
an arbitrary radius R 

(65) 

which we will use to examine the properties of the energy 
functional under shift transformations. 

For all results shown below, the SR calculations used 
as a starting point were performed with an approximate 
particle-number projection before variation within the 
Lipkin-Nogami approach to ensure that pairing correla- 
tions are present in all SR states. Otherwise, pairing cor- 
relations would collapse in the SR state whenever there 
is a large gap in the single-particle spectrum around the 
Fermi surface. 

The dependence of various quantities on axial 
quadrupole deformation is shown in function of the di- 
mensionless deformation of the mass density distribution 
fi2 defined as 

where R = 1.2 A^/'^ fm. 



C. ISO 

As a first example we discuss "'^^O. It has the advantage 
that the density of single-particle levels around the Fermi 
energy is sufficiently low that the impact of the spuri- 
ous contribution brought by each single-particle level to 
the projected energy can be studied separately without 
having them interfere too much. The integration radius 
Rq = 1 \s used until we come to discussing shift invari- 
ance. 



1. Convergence of Operator Matrix Elements 

Before we enter the discussion of the energy functional, 
we demonstrate the convergence of the particle-number 
projection method for observables that are calculated as 
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FIG. 3: (Color online) Dispersion of the proton and neutron 
number of the unprojected SR state and the particle-number 
projected SR using 3, 5 or 7 discretization points of the gauge- 
space integrals as a function of their deformation. For 5 points 
the projected state is sufficiently converged, for 7 and more 
points (not shown) the dispersion cannot be distinguished 
from numerical noise. 



expectation values of the corresponding operators in the 
projected states. In the context of particle-number pro- 
jection, the most sensitive observable is the dispersion 
of particle number (AN^) = (iV^) - (Ar)^, a two-body 
operator that provides a measure for the quality of the 
particle-number projected state as it has to be zero for an 
eigenstate of the particle-number operator. For an (un- 
projected) SR state, (AA^^Ws proportional to its spread 
in particle- number space [70]. One can see in Fig. [3] 
that the Fomenko discretization converges quickly, al- 
ready L ~ 5 gives excellent results for ^^O, and for i > 7 
the dispersion of particle number cannot be distinguished 
from numerical noise. 




FIG. 4: (Color online) Spectrum of poles 



j./v^\ for 



protons (top panel) and neutrons (middle panel), which for 
levels in the vicinity of the Fermi energy resembles a stretched 
and slightly distorted Nilsson diagram. The dashed red line 
at 2 = 1 denotes the radius of the standard integration con- 
tour R — 1. The bottom panel shows the particle- number 
projected quadrupole deformation energy for L — 5 and 199 
discretization points for the integral in gauge space. The in- 
sert shows a close-up of the steps at small deformation. 



2. Regularized PNR Energy 

Unlike any operator expectation value, particle- 
number restored energies do not converge when increas- 
ing the number of discretization points in the gauge-space 
integrals, as already demonstrated in Fig. [1] for the pa- 
rameterization SLy4. Figure |4] shows the projected de- 
formation energy curve of ^^O, now calculated with SIII. 
What appears to be a smooth deformation energy curve 
when calculating it with L = 5, develops steps and dis- 
continuities when increasing the number of discretization 
points to 199, i.e. when one starts to resolve the poles at 
finite close to the integration contour j25j . For exam- 
ple, at small prolate and oblate deformation /32 ~ ±0.15, 
the energy jumps from a lower deformation curve around 
the spherical point to a higher-lying one at larger defor- 
mation. Using a small number of discretization points 
provides a curve that smoothly interpolates between the 



two energy curves distinguished with L = 199. Figured] 
also displays, as a function of the deformation, the poles 
l-^/tl ~ l'^A'/''^A'l that enter uncorrected energy kernels 
for protons and neutrons. We follow Dobaczewski et al. 
[2^ and plot instead of a Nilsson diagram of single- 
particle energies, as divergences and steps appear where 
poles cross the integration contour. Note again that the 
radius of the latter can be chosen to be different from the 
standard value Rq = 1 that is equivalent to the Fermi en- 
ergy 

In Fig. m however, we do not yet make use of the 
freedom to modify the integration contour and use the 
standard values Rp = i?„ = 1. It can be seen that the two 
steps developing at (32 ~ ±0.15 coincide with a pair of 
neutron levels originating from the spherical u dc^ /2+ shell 
that enters the integration contour either at the prolate 
or the oblate deformation. It is noteworthy that the steps 
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FIG. 5: (Color online) Correction for neutrons (top panel) and 
protons (middle panel) and energy gain from projection with- 
out and with correction for '^^O as a function of quadrupole 
deformation for 5 and 199 discretization points for the inte- 
grals in gauge space. The corrected energy gain in indepen- 
dent on the discretization of the integrals when 5 or more 
angles are used. All panels share the same energy scale. 



are not completely sharp even when using L = 199 points 
for the calculation, as can be seen from the markers in 
the insert in the lowest panel. There also is a step at 
/32 = —0.5 that coincides with a pair of proton levels 
from the tt pi/2- shell leaving the integration contour. 
A particular case is the discontinuity at (32 — 0-7 that 
coincides with the crossing of two different pairs of proton 
levels right on the integration contour. 

It is worth noting that no divergence is seen in the PNR 
energy surface displayed in Fig. HI This is at variance to 
Fig. [TJ Indeed, SIII corresponds to a specific functional 
form such that poles at z = zj^ are simple poles. This is 
due to the fact that the trilinear terms entering SIII do 
not contain products of three density matrices referring 
to the same isospin. As explained in Paper III, this leads 
to a finite Cauchy principal value as the poles cross the 
integration circle. Divergences appear only for poles of 
higher order. 

The effects of particle-number restoration on the en- 
ergy is partly masked in Fig. H] by the genuine evolution 
of the energy with deformation. To obtain a clearer pic- 
ture, we show in the lower panel of Fig. [5] the energy gain 
from particle number restoration, obtained as the differ- 
ence between the MR and SR energy functionals for a 



given deformation of the SR state. For a cleaner com- 
parison, the LN correction is removed from the SR en- 
ergy. The steps and discontinuities already seen in Fig. [4] 
appear when increasing L from 5 to 199. The two upper 
panels show the correction 8qq, Eq. ((29)) . separately for 
protons and neutrons. The lower panel also shows the en- 
ergy gain for the regularized PNR energy surface £^^q 
obtained by subtracting the neutron and proton correc- 
tions Sqq from the uncorrected PNR energy £^ for a 
given value of L. The correction has many interesting 
and appealing features 

• The regularized PNR energy E^^g independent 
on the discretization of the integral; it is identical, 
within the numerical accuracy, for i — 5 and 199. 
As a result, only one curve is shown in Fig. [5] 

• The previous result confirms that the entire depen- 
dence of the (uncorrected) PNR energy on the dis- 
cretization of the gauge space integral is contained 
in f^g. 

• Looking separately at protons and neutrons, the 
corresponding correction is largest when a pole 
of a given nucleon species is close to the integration 
contour (i? = 1 here). However, the correction is 
different from zero for the deformations in between; 
i.e. the spurious nature of the poles is also felt when 
being away from divergences and steps. 

• All terms in the energy functional (central, spin- 
orbit, pairing. Coulomb, etc) contribute to f^g, 
with slightly different magnitudes and different 
signs, so one has to strictly correct for all of them. 
This is not unexpected as the source of the spurios- 
ity we focus on here is the weight the matrix ele- 
ments and u"'* are multiplied with in Eq. (|29p . 
not the matrix elements themselves. 

• The correction depends strongly on the deforma- 
tion and will have a non-negligeable impact on the 
topology of the deformation energy curve. The 
regularized energy gain from projection is a much 
smoother function of deformation than the un- 
corrected one, meaning that regularized particle- 
number restoration will provide potential energy 
surfaces with less pronounced structures than un- 
corrected PNR. 

• The correction is of the order of 1 MeV. Of 
course it has to be smaller than the energy gain 
from particle number restoration, which is a few 
McV. For -'^^O however (and when calculated with 
SIII), the spurious contribution to the uncorrected 
energy can be as large as 50 % of the total energy 
gain at some deformations. Also, one MeV error on 
the mass is larger than the targeted accuracy from 
EDF methods. In addition, and as exemplified be- 
low, the correction to the mass varies from nucleus 
to nucleus. 
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FIG. 6: (Color online) Corrected (solid line) and uncor- 
rected (dotted and dashed lines) particle-number projected 
quadrupole deformation energy for ^*0, calculated with L = 5 
and 199 discretization points of the integral in gauge space. 
The corrected curves are identical. 



• The regularized PNR energy gain can be both 
larger and smaller than the uncorrected one. In 
all cases we have looked at so far, however, an in- 
crease obtained from the correction rests always 
very small, while a reduction from correction might 
be quite substantial, but this might not always be 
the case. 

The corrected deformation energy surface of ^^O is 
shown in Fig. [6] together with the uncorrected ones ob- 
tained with L ~ 5 and 199 as was already displayed in 
FigUl It is striking to see that the corrected PNR energy 
surface has less structure than the uncorrected ones; its 
curvature changes now monotonically and the shoulder 
at /?2 = 0.6, that always appears as a secondary mini- 
mum in SR calculations without pairing for oxygen iso- 
topes, disappears completely. The latter does not mean a 
priori that a regularized PNR plus configuration mixing 
calculation will not give a collective state located at this 
deformation anymore as it was obtained for ^^O [7l| and 
^"^O [73 using SLy4. This question needs to be addressed 
in the near future by performing regularized MR calcula- 
tions including quadrupole shape configuration mixing. 



D. Detailed analysis of spurious contributions 

1. Contributions of individual poles 

After discussing the behavior of the contaminated and 
regularized PNR energies of a nucleus as a function of 
its quadrupole deformation, it is instructive to investi- 
gate the contribution of each canonical pair (/x, fl) to 
the unphysical energy £qq that contaminate uncorrected 
MR energies £^ . Formally, each pair of single-particle 
levels provides a spurious contribution through the 



pole at z 

ated with the unphysical poles at finite = ±i|m^/z;^|, 
if the latter are located inside of the integration contour 
of radius Rq. In the end, one can rewrite Eq. (j43p as 
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The total contribution is calculated numerically 
through Eq. ((29| and might depend on the number of 
discretization points L used for the gauge-space inte- 
gral. The partial contribution can be evaluated us- 
ing the analytical expression for the residue of the poles, 
Eq. (|44|) . which docs not depend on the discretization of 
the gauge-space integrals. Finally, £° is equal to when 
1 2^ I > i?, while for |z^| < i? it can be estimated through 
= — As is calculated analytically while £^ 
is obtained numerically, the values obtained for e° might 
not be very precise when \z^ \ ~ R. 

It turns out that only a few pairs of levels located close 
to the Fermi level give a non-zero contribution to £^q. 
The relative size and behavior of these contributions as 
the spectrum of poles changes can be understood by an- 
alyzing Eqns. (|43|) and (|44|) for a few idealized cases. For 
this discussion, the combination of matrix elements enter- 
ing the expression of £^q can be ignored. The values of 
the matrix elements depend of course on the actual pair 
of conjugated states they refer to and thereby scale the 
contribution of a given level to £^q. However, the matrix 
elements do not show a particular dependence on ^ that 
determines the limit of for completely occupied or un- 
occupied levels. Therefore it is sufficient to concentrate 
on the occupation-number dependent weight-factors in 
Eqns. (gS]) and (|44l) . 

Figure [7] separates the various contributions to £qq 
for the three pairs of canonical orbits that originate from 
the spherical neutron ^5/2+ level in ^^O. The top panel of 
Fig. [7] displays the location of the three poles of interest 
on the imaginary axis. Those three pairs of poles are ex- 
plicitly labeled by the jz quantum number denoting the 
projection of the angular momentum on the symmetry 
axis. Other poles are left unmarked. The three other 
panels show e^, and for the three pairs of ^5/2+ 
levels only, as these entirely determine the neutron con- 
tribution to £^Q for the deformations shown^. 

The second panel from the top shows e^. Solid lines 
denote when the corresponding pole is inside the inte- 
gration contour (i?„ = 1 here), while dotted lines denote 
when the pole is outside. Only the former of the two 
contributes to £cq- As ej^ is usually finite when the 



® At large oblate and prolate deformation, the of the other lev- 
els approaching z = 1 are of the same order as those shown, but 
make the plot difficult to read and do not add crucial informa- 
tion. 
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corresponding pole crosses the integration contour, its 
size determines the step left in the PNR deformation en- 
ergy curve. To understand how changes as a function 
of the location of the corresponding pole within the 
spectrum of the other poles, Eq. (|44|l has to be analyzed 
further. The product over v ^ j.i in this expression can 
be estimated by first considering that there arc fc^ pairs 
of levels with \z^\ <^ \z^\, such that their contribution 
to the product can be approximated by 
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For a small number kf of pairs of levels, \z^\ is of the 
same order as \z^\, such that the full factor in the product 
has to be kept. Finally, all remaining levels are such that 
\z^\ <^ \z^\ and the product can again be simplified 
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In practical calculations one works with a limited number 
of pairs of levels kt in the basis. This cutoff, however, has 
no consequence for the contribution from a pair of 
levels {fi, p) below the cutoff, as for all reasonable cutoffs 
the discarded pairs of levels contribute a factor 1 to the 
product in Eq. (|44p . Altogether one obtains 
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where wc assume even particle number A^. Equation (j70p 
allows for the complete explanation of the global behavior 
of £ J seen in Fig. [T] 

First, for a bilinear functional as discussed here, 
is zero whenever the pair of levels (/x, /2) is degenerate 
with another pair (i^jD), i.e. \z^\ = \z^\, as in this case 
the middle product in Eq. (j70|l contains a factor zero. 
In fact, this is a direct consequence of the disappear- 
ance of the pole at z^ in the PNR energy kernel, as the 
dangerous remaining denominator is now canceled by an 
additional factor in the norm kernel^. This alone already 
indicates that the contribution of a given pair of lev- 
els might fluctuate rapidly when the spectrum of poles 




FIG. 7: (Color online) Spurious energy from the single- 
particle orbits that correspond to the spherical neutron 1^5/2+ 
level in ^**0 as a function of quadrupole deformation (see 
text). 



\z^\ changes as a function of a collective coordinate. The 
(— )'^'^ factor in Eq. (j70|l . whose sign depends on the num- 



^ This results holds for any bilinear functional in the density ma- 
trix of a given isospin, even if it is multiplied with the densities 
of the other one. When allowing for higher-order functionals, 
however, a term of order n in the density matrix can generate a 
pole at of order (at most) (n — 1). In order for to be 0, 



one needs the pole at Z/j^ to disappear altogether, which requires 
(n — 1) additional factors from the norm kernel to cancel the 
denominator (ji^ -|- v'^ z'^)~^"~^h Thus, the pair of interest (^, 
p.) needs to be degenerated (at least) with (n — 1) other pairs for 
to be 0. As a consequence, will not be at a simple level 
crossing when working with a trilinear (or higher-order) energy 
functional in the same isospin. 
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ber of pairs of levels kr located below the pair (/x, p,), 
makes to change sign through a crossing with an- 
other pair. Figure [7] contains several such examples. The 
downsloping jz = 1/2"'' substate from the d^^2 spherical 
shell crosses with an upsloping level at large prolate de- 
formation. There, changes its sign as kr changes by 
1 through the crossing. At spherical deformation, where 
the three pairs of ^5/2 levels are degenerate, each of them 
crosses with the two others and kr changes either by 2 (for 
the jz = 1/2+ and jz = 5/2+) or (for the jz = 3/2+), 
such that the corresponding e+ do not change their sign. 
A very particular case is the subsequent crossing of the 
upsloping jz =5/2+ level with two other levels within a 
very small interval around /32 ~ 0.63. As the three levels 
do not cross at exactly the same deformation, e+ changes 
its sign twice in a tiny deformation interval, oscillating 
between values far outside the vertical energy interval 
shown, that cannot be resolved by what appears as a 
single vertical (red) dotted line in the plot at (32 = 0.67. 

Second, let us consider the case of a pair (//, p,) that 
is well separated from all others. Thus, there remains 
only two categories of "other" states in Eq. ([TO)) , k,. pairs 
of levels (C,f) with |z+| ^ |z+| and kt — /ov — 1 pairs of 
levels (A, A) with \z^\ 3> |z+|. One has still to distinguish 
between the two cases where |z+ 1 is larger or smaller than 
1. 

We start with the case jz+j = |M^/iy| > 1 for which 
the factor in Eq. (|70p rapidly converges towards 1 as 
|z+| increases. In this case, the number of pairs below 
the pair (/i, /2) is larger than half the particle number; i.e. 
kr > N/2. For kr = N/2+1, |e+| grows linearly with |z+| 
for |z+| > 1, for kr = N/2 + 2 it grows quadratically etc, 
but always only until it approaches another level, where 
|e+ 1 goes back to as a consequence of the degeneracy as 
described above. After the crossing, however, |e+| grows 
again, although one of the u| w 1 factors in Eq. ([70]) has 
changed into a i;| ^ 1 factor at the crossing. At the 
same time, the number of pairs kr below the pair (/i, fl) 
has grown by one, such that after the crossing there is an 
additional jz+p = m^/w^ factor, that overcompensates 
the effect of the occupation factor from the level just 
crossed, as w| > and v'^ < 1/2 give f| m^/w^ > 1. For 
the simultaneous crossing with more than one level, the 
net effect is the product of the change brought by each 
crossed level. For poles far from the Fermi level, the val- 
ues of e+ can be very large. For example, the e+ of the 
jz ~ 5/2+ level reaches about 550 MeV around (32 = 0.42 
where the corresponding pole |z+| is well isolated in the 
spectrum, drops below zero and rises immediately back 
when it crosses a pair from a higher-lying spherical j 
shell, and quickly rises to values larger than 10^ MeV, 
dropping back to zero right away as the pole crosses the 
next pair, and quickly gaining a value again several or- 
ders of magnitude larger. The sheer size of these values 
that quickly grow beyond any physical scale that appears 
in the EDF description of nuclei clearly shows that e+ 
alone cannot be a meaningful quantity in a well-defined 



theory. The only reason why the e+ of these high-lying 

levels with | z+ 1 3> 1 do not make Sqg incommensurably 
large is that the corresponding poles are outside of the 
standard integration circle and thus do not contribute. 
We will come back to this when discussing PNR with 
shifted contour integrals below. 

For a sufficiently isolated level below the Fermi level, 
|z+| = \Ufj^/Vfj,\ < 1, |e+| also tentatively grows when |z+| 
goes towards 0. This is now a consequence of the fact that 
kr < N/2, such that e+ scales with powers of the inverse 
of |z + 1. At each crossing with a lower lying pair of levels, 
the additional u\ 1 factor is overcompensated by the 
additional |z+|^^ factor from the decreasing number of 
pairs kr below . Again, e+ goes to at level crossings 
and changes its sign depending on the number of pairs 
crossed. 

An important consequence of Eq. (j70p and the discus- 
sion above is that the e+ of an isolated pair is smallest 
when there arc exactly kr = N/2 pairs of other levels 
below it, which is usually the case for a level with its 
pole z+ close to the Fermi level. A side effect is that the 
spurious step due to a pair crossing the integration con- 
tour remains rather small when the latter is chosen as the 
unit circle. This is to put in perspective with the rather 
small spurious steps observed in Fig. [5] and contaminat- 
ing the unregularized PNR energy computed using a unit 
integration circle. We will see in the following that the 
situation would have been more dramatic if we had used 
different contours. 

As discussed in Sec. |Vl poles at finite z+ entering or 
leaving the integration contour are the origin of the spu- 
rious steps in PNR energy surfaces, as the corresponding 
(usually finite) e+ is suddenly added to or removed from 
ScG^ respectively. In the second panel of Fig. [71 contri- 
butions from poles inside or outside the standard integra- 
tion contour of radius R = 1 are plotted as solid or dotted 
lines, respectively, to make this distinction. The third 
panel from the top also shows £+ with solid lines, but 
now only when it actually contributes to ScG- dot- 
ted lines represent — such that the distance between 
the curves for e+ and — e° provides the total contribution 

from the pair (/-t,/i) to £qq ^■ 

The results for the neutron levels depicted in Fig. [7| 
suggest that £+ converges towards — when z+ goes 
to zero, i.e. for deeply-bound levels far below the Fermi 
energy, such that the total contribution is zero for 
deeply-bound levels. When z+ approaches the Fermi en- 
ergy from below, e+ and — eJJ slowly grow apart. Still, for 
all examples we have looked at, e|J and e+ remain of simi- 
lar size, but opposite sign, and have a similar dependence 
on deformation around the Fermi energy, z+ « 1. They 
do not cancel exactly when the pole at z+ approaches 



The spikes of at the deformations where the contribution from 
to S^Q jumps to are of numerical origin. 
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the Fermi level but the difference between and — e°, 
remains much smaller than the individual contributions 
and provides the finite and smoothly varying spurious 
energy Eqq between the steps. For levels far above the 
Fermi level, e° goes to zero. Also, the pole is beyond 
the integration contour and does not contribute to 
either. Consequently, levels far above the Fermi en- 
ergy do not contribute to £qq for standard integration 
contours at i?g = 1. 

The behaviors described above can be understood as 
limiting cases of the factor u^w^ times the contour inte- 
gral in Eq. (j43|) . Omitting unimportant prefactors, one 



obtains for \z^\ 



0, that is for 



and vf, w 1, that 



Indeed, right at the crossing behaves as 



= K I't (4-4 [m, C] - 2 c%_,[fi, C] + c% [m, C]) (75) 

where c^[/i, ^] denotes a modified norm obtained by re- 
moving the contributions of both pairs (/i, p.) and (C, (^) 
from the usual norm kernel 
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while for |z^| — * oo, that is for ~ 1 and « 0, one 
has 
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where the cat denote in both cases the amplitudes of the 
normalized projected states with particle number N in 
the SR state, Eq. ([42]) . all of which are usually non-zero 
and independent of fi. The key element to obtain both 
limits is that the integral over the gauge angle becomes 
simply proportional to « 1 or « 1, respectively. 
As a result, the prefactor f ^ dominates and drives 
towards zero in both cases. As a consequence, one indeed 
finds as a general rule that 
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as suggested by the numerical results in Fig. [T] 

Unlike e^, the contribution £° to the physical pole at 
z = is not a priori suppressed for degenerate levels and 
might have a non-zero value. For deep-hole states, this 
seems contradictory with the previous proof that = 
goes to zero. In fact, when the pair (/i, /2) crosses 
another one (CjC)i not only the pole at is removed 
but the residue of the pole at z = is strongly affected 
by the disappearance of the corresponding denominator. 
As a result, £° also goes towards zero as goes to zero. 



Considering either rather deep-hole or highly-lying 
single-particle states, the prefactor (m^ v^)'^ appearing in 
Eq. (ffS]) makes = e° to be small. 

The bottom panel of Fig. [7] shows the total contri- 
bution Efj. of each selected pairs to Sqq- One can now 
clearly see that there is more to the spurious energy than 
just the steps and the divergences (the latter of which do 
not appear for the particular functional used here). The 
poles z^ associated to the jz = 3/2^ pair remain out- 
side the integration contour for all deformations. Thus, 
it does not produce a step as the corresponding never 
contributes to £qq. Still, this level gives a small con- 
tribution to the spurious energy through the pole at 
z = 0, which happens to be slightly larger for prolate 
deformations than for oblate ones. 

Starting on the oblate side, only the pole at z = con- 
tributes at first to the spurious energy from the jz = 1/2"*" 
pair of levels. The corresponding ej^ increases slowly from 
zero with increasing (32- The moment the corresponding 
poles z^ enter the integration contour at /32 = 0.15, 
suddenly contributes to the spurious energy. We already 
saw that the finite value of at this point determines 
the step. As approaches — e° when the 1/2+ levels 
become more and more occupied with increasing prolate 
deformation, the total contribution of the 1/2+ pair to 
now decreases after the step. With the total contri- 
bution increasing on one side of the step and suddenly 
decreasing on the other, the curvature of the spurious 
energy is different on both sides of a step. As a conse- 
quence, removing £^q modifies the improper curvature of 
the uncorrected deformation energy surface that is vis- 
ible in Fig. [6l even when the steps themselves are not 
numerically resolved. The regularized deformation en- 
ergy surfaces show much less structure; in ^*0 to the ex- 
treme that the curvature of the corrected energy surface 
is now positive for all deformations as shown in Fig. [6l 
The contribution from the 5/2"*" levels to the spurious 
energy behaves very much as the one from the 1/2+ lev- 
els with oblate and prolate sides exchanged. The sum 
of the three individual contributions gives the neutron 
correction shown in the top panel of Fig. [5] for L ~ 199. 

We have seen that for a bilinear functional, the steps 
are always the consequence of a pair of single poles 
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FIG. 8: (Color online) Same as Fig. [T] but for the proton 
1/2^ and l/2~ levels that give the dominant contributions to 
the proton correction at prolate deformation and cross at the 
Fermi energy at /32 = 0.67. 



crossing the integration contour and have the size of the 
corresponding at that crossing. The steps cannot add 
up for a biUnear functional as, for degenerate poles with 



zp directly lead to = ejj, = 0. This 



does not mean, however, that there is no spurious con- 
tribution to the PNR energy when two poles cross the 
integration contour simultaneously, as the corresponding 
ej^ and e^, are in general non-zero. In fact, Fig.[5]demon- 
stratcs clearly that, in our calculation of with SIII, 
the spurious energy £qq is largest exactly where two pro- 
ton levels cross at the Fermi energy at /?2 = 0.67. The 
contribution of these two pairs of levels to S^G' which 
also happen to be the only proton levels that have a 
finite contribution at prolate deformation, is analyzed 
in Fig. [8l There is a number of interesting differences 
with Fig. [T] (i) The contribution does not vanish at 
spherical shape for the l/2~ levels for a bilinear func- 
tional. Indeed, the spherical pi/2 shell is only doubly 
degenerate, which does not suppress the corresponding 
e^. In fact, only si/2 andpi/2 levels with poles below 
the integration contour provide non-zero at spherical 



shape, (ii) Both pairs cross right at the Fermi energy at 
j32 = 0.67. For the standard choice = 1, their poles 

thus cross on the integration contour. As a result, 
from both pairs are zero, and change sign at the crossing, 
(iii) The derivative of is not zero for both pairs when 
they simultaneously cross the Fermi energy. By contrast, 
ej^ slowly approaches zero such that the total contribu- 
tion is quite large for the two proton pairs. When the 
poles z^ approach the integration contour from below, 
the distance between and e|J grows for both pairs. 
After the poles have crossed the contour, only the £° 
contribute. Finally, the total contribution from each 
pair that crosses with another at the integration contour 
is largest at the crossing, and decreases towards zero on 
both sides. The sum of the two individual contributions 
gives the proton correction shown in the middle panel of 
Fig.Ofor L — 199; all other proton levels are too far away 
from the Fermi level to provide any visible contribution. 

One can take advantage of the fact that only a very 
limited number of levels actually contributes to £qq in 
order to reduce the numerical effort. Evaluating the nec- 
essary matrix elements v'''^ and v'^'^ only for those levels 
for which the weight is significantly different from zero is 
particularly welcome for the expensive contribution from 
the Coulomb interaction. 



2. Shift Invariance 

In their recent paper, Dobaczewski et al. [2^ pointed 
out that the (uncorrected) PNR energy density func- 
tional is not shift invariant, i.e. PNR energies depend 
on the radius chosen for the contour integral in the com- 
plex plane. As outlined in Sees. IV D] and fVEl the source 
of violation of the shift invariance is the contribution 
from the poles at z^ inside the integration contour Cn 
to the spurious energy S^q in Eq. . Each time a pole 
z^ enters or leaves the integration contour when chang- 
ing its radius, Sqq changes by the amount e^. This is 
illustrated in Fig. [9] for ^^O at /32 = 0.371. The radius 
of the contour used for neutrons is held fixed at i?„ = 1, 
while the radius of the contour used for the protons is 
varied. The three steps visible in Fig. ^ correspond to 
the three proton poles located at 0.1 < < 10 visible 
in Fig. m for the deformation of interest. 

An interesting feature of the steps is that their size 
grows as the integration contour is shifted away from 
Rp = 1 [1^, i.e. away from the Fermi level. The reason 
is easy to understand from the discussion of Eq. (jTO]) 
given in the previous section: increases as \z^\ moves 
away from 1 (as long as it is separated from other poles) 
and as the difference between the number kr of pairs of 
states below (/x,/z) and half the number of particles N/2 
one is restoring grows. 

Using the small number of L = 5 discretization points 
for the gauge-space integral does not resolve the steps 
in the uncorrected PNR energy; only with much larger 
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FIG. 9: (Color online) Projected energy for at the de- 
formation /32 = 0.371 as a function of the radius Rp of the 
integration contour calculated without and with correction 
using 5 and 199 angles. The energy scale on the left gives the 
absolute energy, the scale on the right the energy gain from 
projection. The insert magnifies the curves around -Rp = 1. 
The regularized PNR energy in independent on the discretiza- 
tion of the integrals when 5 or more angles are used. The in- 
tegration contour for projection on neutron number is _R„ = 1 
in all cases. 



L one obtains sharp steps. By contrast, and as seen in 
Fig. [21 the regularized PNR energy is constant within a 
numerical precision of the order 1 keV as Rp is modified 
and L increased beyond 5. 
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FIG. 10: (Color online) Weight c%{Rp = 1) = 
of the normalized proton-number projected states in the 
SR HFB state (upper panel), the weighted spurious energy 
Cz{Rp = ^) ScaiRp — 1) (middle panel), the non-regularized 
weighted PNR energies c%{Rp = 1) £^ {Rp — 1) and regular- 
ized c%{Rp = 1) £b.eg{Rp — 1) (lower panel). All results are 
obtained using the same SR state calculated for at a 
deformation of P2 ~ 0.371 as auxiliary state. The neutron 
number is not restored. 



3. Distribution of weighted PNR energies 

As a next step, we analyze how the spurious energy 
£qq{R) affects the distribution of non- normalized PNR 
energies cj^{R) E'^ (R) and c^(l)£^(i?) as a function of 
the particle number one restores. Of course, restoring 
other particle numbers than the one that the underly- 
ing SR state was constrained in average to is not very 
useful for practical applications. The purpose of the ex- 
ercise, however, is to shed further light on the nature of 
the spurious energy £^q[R), especially through testing 
sum rules associated with such a decomposition over N. 
For the latter test to be meaningful, and as motivated 
in Sec. IV El it is essential to include zero and negative 
particle numbers in the analysis. 

Starting with a SR calculation for ^'^O, the average 
proton and neutron number are small enough that non- 
zero values of the quantities of interest for negative par- 
ticle numbers can be unambiguously detected in the tail 
of the distribution when performing a numerical calcu- 
lation. Of course, a SR state with even number-parity 
quantum number, as assumed here, can only be projected 
on even particle number, such that the weight c^{R) and 
any operator matrix elements are obviously zero for odd 
N. In addition, the contributions to f ^(i?) from the spu- 



rious poles, see Eq. ([33|), and from the physical pole^ are 
zero for odd N when restoring particle number from a 
SR state with an even-number parity quantum number. 
As a consequence, we can limit ourselves here to looking 
at even particle numbers. 

For the sake of transparency, and to avoid double sums 
over TV and Z as well as the interference of the corre- 
sponding terms when analyzing the sum rules, we limit 
ourselves to the restoration of proton number in this sec- 
tion and in the following one. We start with the same SR 
state calculated for ^^O with f32 = 0.371 as in Fig. [5] but 
without restoring neutron number, which is constrained 
to an average value of iV = 10. The restoration of proton 
number is performed using L ~ 199 integration points. In 
what follows, we discuss the interaction part of the EDF 
only, i.e. the EDF without kinetic energy and without 
the one-body center of mass correction used in connec- 
tion with SIII. Both are expectation values of one-body 
operators and therefore free of spurious contributions. As 



The Laurent series centered at 2: = of the integrand in Eq. II35I I 
does only contain even powers for odd N. As a result, such a 
pole does not contribute to £^ (R). 
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before, the Coulomb exchange term is omitted from the 
energy functional. 

First, we discuss the standard case with an inte- 



gration contour at Rp 



The upper panel of 



Fig. [TO] displays the distribution of the weights c| {Rp ~ 
1) = Eq- dm), of the normalized proton- 

number projected states in the SR state. As expected, 
c^{l) is peaked at Z = 8 and falls off quickly to numeri- 
cal noise. Components with Z > 14 and Z < 2 cannot be 
numerically distinguished from zero. In the former case 
and for Z = it is a consequence of these proton num- 
bers being too far from the average proton number such 
that c^il) becomes too small to be distinguished from 
zero within the numerical precision of our code, while 
for Z < the proton-number projected states l^"^) are 
strictly zero for analytical reasons. 

The lower panel of Fig. [TO] shows the interaction part 
of weighted PNR energies before and after applying the 
regularization method. The distribution of absolute val- 
ues of c^(l) f ■^(1) does not follow the distribution of the 
weights c\(l) displayed in the upper panel. Instead, it 
has a long tail that spreads visibly to Z = — 20 and 
Z = 34, before it cannot be distinguished from numerical 
noise anymore. In these tails, the values of c\{\) £^ {\) 
have alternating signs, which is clearly unphysical. Only 
the regularized quantity c\{l) £^^q{1) falls off in the 
same manner as c\{l) does and is numerically zero for 
Z < 0. This underlines again the spurious nature of £qq, 
that is shown separately in the middle panel of Fig. [TO] In 
the present example, c\{l) £qq{1) has alternating signs 
throughout the entire interval of Z. This must not al- 
ways be the case; in some other examples we have looked 
at, this happens only for particle numbers that are a at 
least a few units away from the average particle number 
of the underlying SR state. 

For Z < 0, non-zero c^(l)£'^(l) are entirely spurious 



with £^{1) 



1); i.e. they originate entirely from 



spurious poles at finite z*. The same situation applies 
to the tail of the distribution of c|(l)£^(l) for large 
positive Z. 

As a second example, we show in the three upper 
panels of Fig. [TT] the same quantities as in Fig. [TOl 
but obtained employing an integration contour of radius 
Rp = 2.5. By contrast to before (Rp = 1), the poles 
from the 1/2"*" substate of the tt 1^5/2+ shell are located 
inside the integration contour, see Fig. [5] As a result, 
the spurious contribution from those poles increases 
£^ by about 4 MeV when projecting on Z = 8 using a 
non-regularized functional, see Fig. [9] We analyze now 
if and how the energy restored on other proton numbers 
are affected compared to using Rp = 1. 

Compared to Fig. [TO] the distribution of weights 
c|(2.5) is distorted by the additional Rp = (2.5)^ fac- 
tor such that absolute values change by several orders 
of magnitude, and the maximum of the distribution is 
shifted to Z = 10. The main difference to the case us- 
ing the standard integration contour Rp = 1 is that the 
distribution of the spurious energy c^(2.5) £qq{2.5) is 
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FIG. 11: (Color online) Weight c%{Rp = 2.5) = \$R^)f 
of the normahzed proton-number projected states into the ra- 
dially shifted SR HFB state at Rp = 2.5 (upper panel), the 
weighted spurious energy c%{Rp = 2.5) ScciRp (upper 
middle panel), the non-regularized £^ {Rp = 2.5) and regular- 
ized f^B<3(i?p = 2.5) PNR energies weighted by c%{Rp = 2.5) 
(lower middle panel) and by c%{Rp = 1) (lower panel). All 
results are obtained using the same SR state calculated for 
at a deformation of [32 = 0.371 as auxiliary state. The 
neutron number is not restored. 



distorted in a different manner than the distribution of 
the norm, such that it falls off quicker for Z > 8, but 
much slower for Z < 8, including negative Z. Again, 
only the distribution of the regularized MR energy func- 
tional £^^q{2.5) follows that of the weights c^(2.5). 

The lowest panel of Fig. [Tl] shows the contributions 
to the non-radius-weighted energy sum rule discussed 
in Sec. IVE2I The distribution c%{l)£^{Rp) is even 
more distorted than for the contributions to the radius- 
weighted sum rule shown in the panel above. For Rp > 1, 
4(1) £^{Rp) falls off much quicker than c|(i?p) £^{Rp) 
for Z > 8, but much slower for Z < 8. For negative val- 
ues of Z the missing factor (2.5)^ makes c|(l)£^(i?p) 
to grow so fast that it will be impossible to safely evalu- 
ate numerically the sum rules including negative particle 
numbers. 
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For Rp < 1, a case not discuss here, the situation 
is reversed such that c\{\)8^{Rp) falls off faster than 
c%{Rp)£^{Rp) for Z < 8, but slower for Z > 8, now 
with the impossibility to safely evaluate the sum rule 
when including positive Z . 

To summarize, the contamination of the PNR EDF 
by spurious contributions originating from the use of 
the GWT impacts the decomposition of the (shifted) 
SR functional energy (kernel) into weighted PNR en- 
ergies with different particle numbers such that energy 
is shifted out of the physical subspace corresponding to 
positive particle numbers. The impact of this finding on 
the fulfillment of basic sum rules is examined in the next 
section. 



4-. Sum rules 

Now we turn to the sum rules, which arc obtained by 
summing the weighted PNR energies shown in Figs. [TUl 
and 111! Numerical summation is performed on a subset 
of even proton numbers in the interval — 20 < Z < 40. 

Again we begin with the case i?p = 1, for which the 
radius-weighted and non-radius-weightcd sum rules are 
identical. The SR encrgy^° that sets the reference is given 
by 

= -410.3403 MeV. (77) 

In agreement with Eq. ([52|) . the sum of c|(l) 5^(1) over 
positive and negative Z reproduces this value better than 
0.1 keV 

4(1)5^(1) = -410.3403 MeV. (78) 

Z=-oc 

When limiting the sum to "physical" proton numbers 
Z > 0, however, we obtain instead 

J2 4(1) ^^(1) = -410.3550 MeV . (79) 

Z>Q 

With 14.7 keV, the numerical difference between Eq. ([78]) 
and ([75]) . which constitutes the breaking of the physical 
sumrule, is quite small. Using the standard integration 
contour of Rp = 1, we find similar values for other de- 
formations in ^^O, whereas for heavier nuclei this quan- 
tity becomes rapidly suppressed, such that it cannot be 
unambiguously detected in a numerical calculation any- 
more. 

The largest individual sum-rule breaking contribution 
is that for Z = 0, for which wc obtain 

c|(l) £^=°(1) - 4(1) £g^'>{l) = 0.0189 MeV , (80) 



We recall that quoted energies are without the kinetic and centcr- 
of-mass correction energies. 



which is slightly larger than the entire sum over Z < 0. 
This is not unexpected in view of the alternating signs of 
the contributions pointed out in the previous section. 

For Z < 0, non-zero c^(l) are of course entirely 

spurious, such that they equally contribute to the sum 
rule of c\{l) Eqq{\). For Rp — 1, the right-hand-side of 
Eq. ([S3]) is zero, such that the sum of c%{l) £^q{1) over 
all Z is zero as well, which we do find numerically 

-)-oo 

4(1) '^^cg(I) = 0-0000 MeV. (81) 

Z=-oc 

Although the alternating sign of 4(1)^cg(1) with Z 
indicates that a cancelation effect is at play, the result 
of Eq. ([5T|) is not so obvious when looking at the middle 
panel of Fig. [TOl Summing up 4(1)^cg(1) ^'^^ positive 
values of Z gives —0.0146 MeV, which precisely is the 
difference between Eqns. ([77)) and ([75]) . 

The regularized energy 4(1) '^fi_EG(l) numerically 
zero for Z < as any meaningful particle-number re- 
stored observable should be. The same holds for those 
large positive values of Z where > 0. As a conse- 
quence, the sum over 4(1)^K£;g(1) '^^^ be limited to 
"physical" particle numbers. The numerical value for 
this sum 

-t-oo 

E 4(l)f|i=;G(l) = E4(l)f|i.G(l) 

Z=-oo Z>0 

= -410.3403 MeV (82) 

gives back the SR energy. Eq. ([77]) . within 0.1 keV as 
expected from Eq. ([CT|) . 

When shifting one of the states to Rp — 2.5, the norm 
kernel is (<&i|$2.5) = 2816.9760, and the corresponding 
transition energy kernel is £[2.5] = -830.2386 MeV. The 
reference for the radius-weighted sum rule is thus pro- 
vided by 

£:[2.5] ($i|$2.5) = -1266844 MeV, (83) 

where we limit ourselves again to seven digits. Summing 
4(2.5) £■^(2.5) over all Z reproduces this value with the 
same precision, while summing over positive Z only gives 
— 1266546 MeV, which differs from the above value by 
—298 MeV, which is of similar order as in case of the the 
unshifted integration contour. 

In the case of shifted contours, the non-radius- weighted 
sum rule is more interesting to look at. As became clear 
from the bottom panel of Fig. [Tl] the sum over all Z 
cannot be evaluated numerically. Let us anyway focus on 
the sum rule over positive Z only; i.e. Eq. (j62p . In this 
case, summing 4(1) ^^(2.5) gives -309.4217 MeV which 
indeed decomposes into £[p,K,K*] = —410.3403 MeV 
plus the sum-rule breaking term obtained (either nu- 
merically or analytically through Eq. (|63p ) by sum- 
ming c%{l) £qq{2.5) over Z > and which equates 
-1-100.9189 MeV. Thus, one realizes that the most essen- 
tial non-radius-wcightcd sum rule performed over physi- 
cal components (Z > 0) is broken and not shift invariant. 
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FIG. 12: (Color online) Weight of the normalized state pro- 
jected on various values of Z in the SR vacuum (top panel) 
and decomposition of the energy into Z components for three 
different radii of the integration contour for protons (bottom 
panel) for at a deformation of (32 = 0.371. All states 
are projected on the same neutron number TV = 10 with an 
integration contour of radius R„ — 1, using L — 199 integra- 
tion points for both protons and neutrons. Corrected PNR 
energies are the same for all values of Rp within numerical 
accuracy. 

In particular, the breaking term can be very large as soon 
as the integration radius differs from 1. Of course, the 
small (non-zero) value of that sum-rule breaking term 
obtained from using the unit circle as an integration con- 
tour has masked the contamination of energy functionals 
with spurious terms for many years. Indeed, practition- 
ers naturally interpreted that very tiny breaking as due 
to numerical noise. 



5. Energies of physical systems 

After looking into the contributions to the sum rules, 
we now turn our attention to normalized PNR energies 
pertaining to the physical subspace, i.e. addressing only 
those particle numbers that give a non-zero norm. Fig- 
ure [12] shows PNR energies (now again completed by ki- 
netic energy and cm. correction) for three values of the 
integration contour radius Rp. With each step in the un- 
corrected projected energy of the Z = 8 component seen 
for Rp — 1.9 and Rp = 8.2 in Fig. [9l the energy of all 
other Z components changes as well. For each radius of 
the integration contour there is at least one Z compo- 
nent that has an obviously unphysical uncorrected PNR 
energy. 

The breaking of the physical sum rule for the non- 



regularized PNR EDF discussed above is much smaller 
then the energy scale of the changes we observe in Fig. [H] 
when shifting Rp. Still, we can argue with the help of the 
sum rules for the regularized and non-regularized PNR 
EDF that any small spurious energy in a Z component 
with large weight might have to be compensated by a 
very large spurious energy in a Z component with small 
weight, as it happens in Fig. [T^lfor the Z = 12 compo- 
nent at Rp — 1.0 and the Z = 6 component at Rp = 4.0. 
As a consequence, the moderate energy scale found for 
the spurious energy along a deformation energy surface 
when projecting on the same nucleon number that SR 
vacua were constrained to does not apply to the spurious 
energies entering other Z components. While this usu- 
ally has no particular consequences for particle restora- 
tion calculations where one is in most cases interested 
in projecting out the one particle number that the SR 
HFB state was constrained to and which can be expected 
to have a comparatively small contamination of spurious 
energy, the spurious redistribution of energy might seri- 
ously compromise angular-momentum restoration, where 
one is often interested in producing the entire spectrum 
of low-lying states. 



E. ^^Kr 

With the next example ''^Kr, a medium heavy nucleus 
located in a region of shape coexistence, we examine 
how the spurious contributions to the particle-number 
restored energy evolve when increasing the density of 
single-particle levels. This nucleus is one out of the series 
of neutron-deficient Kr isotopes that were recently stud- 
ied [ill with GCM mixing of quadrupole deformed axial 
particle-number and angular-momentum restored states 
using SLyG.^i 

Figure [T3] shows the location of the poles at for 
protons and neutrons, the energy gain from PNR and 
the absolute PNR energy as a function of quadrupole de- 
formation, both with and without correction and both 
calculated with L = 9 and 99 discretization points of 
the gauge-space integrals. We have checked that all ob- 
servables calculated as operator matrix elements are con- 
verged for L = 9. The main difference to ^^O is the 
much larger overall density of poles. This has two con- 
sequences, (i) It increases the number of poles crossing 
the integration contour when deforming the nucleus and 
thus the number of steps, (ii) Poles crossing the Fermi 
level are hardly isolated from other poles which limits the 
size of the steps through the factors entering the middle 



^ The deformation energy surface obtained with SIII also displays 
shape coexistence, although its topography is quite different from 
the one obtained with SLy6. With SIII, the deformed minima 
are much more pronounced and lower in energy compared to 
the spherical one. However, this is irrelevant for the present 
discussion. 
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FIG. 13: (Color online) Spectrum of poles = \u^/vi_t\ 
for protons and neutrons, the uncorrected and corrected en- 
ergy gain from projection and the particle-number projected 
quadrupole deformation energy for L = 9 and 99 discretiza- 
tion points of the integral in gauge space for '^''Kr. 



product in Eq. ([70)) . As a consequence, most of the steps 
visible in Fig. [13] are much smaller than those found for 
^^O in Fig. [H Notable exceptions are the ones on both 
sides of the prolate minimum at P2 ~ 0.43, which in- 
deed correspond to the crossings of proton levels that arc 
well separated from other poles. The correction is not of 
the same magnitude in the various minima; in fact, the 
variation of the correction between the various minima 
is of the same order as the difference in total energy of 
the latter. Correcting for spurious energies might have a 
visible impact on the excitation spectrum of this nucleus 
obtained from a GCM mixing over quadrupole shapes of 
particle number and angular momentum restored states. 
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FIG. 14: (Color online) Spectrum of poles Zfj, — \Ufj_/v^\ for 
protons and neutrons, the correction to the particle num- 
ber restored EDF separately for protons and neutrons, the 
uncorrected and corrected energy gain from projection, and 
the particle-number projected quadrupole deformation energy 
without and with correction for L = 13 and 99 discretization 
points of the integral in gauge space for '^^^Pb. 



'Pb 



As the last example, we present in Fig. [14] results ob- 
tained for -'^^^Pb, a nucleus exhibiting triple shape co- 
existence of spherical, oblate and prolate states studied 
earlier in Refs. [6l,[6l] in a method that includes particle- 
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number restoration using the Skyrme EDF SLyG"'^^. In 
this heavy nucleus, the number of neutron poles in 
the vicinity of the Fermi level is even larger than for 
''^Kr. When crossing the standard integration contour 
Rn = 1, those poles generate many steps which are, 
however, almost always of tiny size due to the close- 
ness of other poles; the sole exception being the step 
at /?2 = 0.4. This is different for protons. As a conse- 
quence of the magic proton number Z — 82, the density 
of proton poles around the Fermi level is quite low for 
most deformations, such that the few proton poles that 
cross in these regions have a much larger impact. This 
is illustrated by the second panel in Fig. [14] that shows 
the correction £qq separately for protons and neutrons. 
The narrow peak at small oblate deformation (32 = 0.11 
is not a divergence, but stems from the crossing of two 
proton levels at the Fermi energy in analogy to the struc- 
ture found in ^^O around /?2 = 0.67. In both cases the 
double-crossing is a direct consequence of the shell clo- 
sure: With all other levels being too far above or below to 
have occupation numbers significantly different from or 
1, the constraint on the average particle- number dictates 
that two pairs of levels in the gap have an occupation of 
vj^ = 1/2 simultaneously. Interestingly, the uncorrected 
deformation energy curve does not change much when in- 
creasing the number of integration points from L = 13 to 
99. As for ^^O and ^^Kr, the correction varies strongly 
with deformation, has a different value in the various 
minima, and, most importantly, is on the same energy 
scale as the energy difference between those minima. 

VII. SUMMARY, CONCLUSIONS AND 
OUTLOOK 

In the present paper, we introduce the notion of spuri- 
ous self-pairing. It appears as a generalization of spurious 
self-interaction processes, a well-known problem in elec- 
tronic density- functional theory 0, [H, [60, [fill, [HI] , to sys- 
tems with pairing correlations that are modeled within 
EDF approaches using independent quasiparticle BCS 
states as auxiliary states of reference. Self-interaction 
and self-pairing processes appear for any energy func- 
tional that uses different vertices in the particle-hole and 
particle-particle channels, and/or not fully antisymmet- 
ric vertices; e.g., as due to density-dependencies. Neither 
self-interaction nor self-pairing appear when the many- 
body energy is strictly calculated as the expectation value 
of a Hamilton operator. Both are a price to pay when 
replacing the exact nuclear many-body problem by a sys- 
tem of coupled one-body problems in a EDF calculation, 
modeling higher-order in-medium correlations through a 



The deformation energy surface obtained with SIII is at variance 
with the experimental finding that the ground state is spherical 
with low-lying prolate and oblate bands seen as excitations [68L 
|69|| . However, this is irrelevant for the present discussion. 



simple energy functional depending on one-body den- 
sities and currents. On the single-reference level, self- 
pairing gives a spurious contribution to the pairing field 
(and therefore influences all quantities it determines) and 
to the total binding energy. 

Energy density functionals extended to perform multi- 
reference calculations, i.e. symmetry restoration or 
GCM-type configuration mixing, also contain unphysi- 
cal contributions: First, the previously discussed self- 
interaction and self-pairing processes that continuously 
extend from SR energy functional to off-diagonal energy 
kernels, as well as a second and much more dangerous 
category of spuriosity that appears when the off-diagonal 
kernels are evaluated on the basis of the generalized Wick 
theorem. The use of a Wick theorem to evaluate a func- 
tional energy kernel that does not originate from a gen- 
uine Hamilton operator is not justified. Relying on the 
generalized Wick theorem to construct off-diagonal func- 
tional energy kernels has the unexpected particularity 
to provide previously discussed self-interaction and self- 
pairing contributions with unphysical weights that con- 
tain poles leading to divergences [3| and steps in the 
energy [2^. The latter have been noticed recently in 
the context of particle-number restoration whenever a 
single-particle level crosses the Fermi energy. As demon- 
strated in Paper I , the weights of self-interaction and 
self-pairing terms obtained on the basis of the standard 
Wick theorem are different and do not present any prob- 
lematic contributions. This feature can be exploited to 
unambiguously isolate the dangerous spuriosities and set- 
up a correction scheme that regularizes unphysical di- 
vergences and steps in MR energy kernels [2^. In the 
present paper, we have applied this correction scheme to 
the simplest and formally most transparent MR case of 
particle-number restoration after variation. 

The complex-plane analysis performed in the present 
work reveals that each conjugated pair of single-particle 
levels (/i, p.) provides an unphysical contribution to the 
physical pole at z = 0, in addition to generating unphys- 
ical poles at = ±i\Ufj,\/\v^\. The latter cause the steps 
as they cross the integration contour [2^. The unphysi- 
cal poles are also at the origin of the breaking of the shift 
invariance of PNR energies [2^ . However, removing only 
the contribution from the poles at z^ to the energy func- 
tional kernel leads to unphysical results. Instead, the 
spurious contribution from a given pair of single-particle 
levels to the pole at z = has to be removed simultane- 
ously, as both are very large, of opposite sign, and nearly 
cancel. 

The correction scheme proposed in Paper I docs indeed 
remove both contributions; thereby it eliminates the di- 
vergences and steps and restores the shift invariance of 
PNR energies £^ as well as standard sum rules that they 
can be expected to fulfill. The correction to £^ is of the 
order of 1 MeV, and in most cases reduces the energy 
gain from PNR. On the one hand, the correction is suffi- 
ciently small that PNR EDF results published earlier are 
not meaningless. On the other hand, in extreme cases the 



32 



correction might be as large as 50 % of the energy gain 
from PNR, which casts some doubt on the rehabiUty of 
pubhshed calculations performed within the EDF frame- 
work. The correction is also of the same order as the rms 
error of the mass residuals reached with the best avail- 
able particle-number restored EDF mass fits [l^. The 
correction varies rapidly with deformation and affect sig- 
nificantly the structure of complex nuclei presenting soft 
deformation energy surfaces and coexisting minima. 

In the present paper, we do not attempt to correct for 
the "true" self-interaction and self-pairing processes that 
contaminate the single-reference energy density function- 
als. This amounts to modify the underlying functional 
which we postpone to later works. In addition, a self- 
consistent correction is very cumbersome, as documented 
in the literature for self-interaction in the context of elec- 
tronic DPT [i, M, S EH, Ei] . 

Particle-number restoration is not the only type of 
MR-EDP calculation where using the GWT as a ba- 
sis to construct non-diagonal functional energy kernels 
causes problems. In fact, any symmetry restoration or 
GCM-typc configuration mixing calculation is expected 
to be contaminated with similar spurious contributions; 
e.g., anomalies were encountered in Ref. [sot in angular- 
momentum restoration calculations of cranked states 
without pairing and using a Skyrme EDF. The correction 
scheme proposed in Paper I can be applied to any type of 
MR-EDF calculation. However, all others but particle- 
number restoration require the numerical construction of 
the canonical basis of the Bogoliubov transformation con- 
necting the two different quasi-particle bases associated 
with the two vacua entering the construction of the func- 
tional energy kernel [1^. Work towards the numerical 
implementation of such a scheme is underway. 

In the present study, we have limited ourselves to par- 
ticle number restoration after variation, where the cor- 
rection can be subtracted from energy kernels a posteri- 
ori. With variation-aftcr-symmetry-rcstoration EDF cal- 
culations becoming available ^26, , 52J , and the variational 
equations sometimes running into the divergences [25j . 
setting up a correction scheme for those variational equa- 
tions becomes an important issue and will be addressed 
in a forthcoming study [8l| . 

The correction proposed in Paper I and discussed 
in the present one is limited to energy functionals de- 
pending on integer powers of the density-matrices. Most 
functionals used in the literature, however, have a den- 
sity dependence of non-integer power, both in the func- 
tional modeling the effective strong interaction and as 
an approximate Coulomb exchange term. Compared to 
the functionals discussed here, such non-integer powers 
of the density matrix pose two additional types of diffi- 
culties when extended to non-diagonal energy kernels on 
the basis of the GWT: (i) as transition densities are com- 
plex, taking their fractional power is ambiguous (25j , and 
(ii) there is no well-defined basis at present to remove the 
spurious branch cuts that are generated by such terms. 
Both points are illustrated and examined further in Pa- 



per III of this series [27[ . 

In our opinion, the particular difficulties of functionals 
with non-integer density dependencies constitute a strong 
motivation to construct energy functionals with integer 
powers of the densities only in view of performing mean- 
ingful MR-EDF calculations in the future. At present, 
there are no such non-relativistic functionals of high per- 
formance. Relativistic functionals have been constructed 
along these lines recently [t^ with a different motiva- 
tion, and have already been used in PNR-EDF calcula- 
tions [20| . The construction of correctable energy func- 
tionals for multi-reference applications becomes an ur- 
gent task for the future. A particular problem will be 
to find a suitable functional for the Coulomb interaction, 
as using the exact exchange term is incommensurately 
expensive in multidimensional MR-EDF calculations. 
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APPENDIX A: THE ENERGY FUNCTIONAL 

The energy is given as the sum of the non-interacting 
kinetic energy, the Skyrme energy functional that mod- 
els the strong particle-hole interaction, a pairing func- 
tional that models the particle-particle interaction and 
the Coulomb energy functional 



So 



ioulomb 



(Al) 



The kinetic energy is the mean value of a one-body op- 
erator; hence it does not pose problems. From the point 
of view of establishing the correction to the MR energy 
kernel, we identify in the following 



£PP 

£PPP 



cPP I cdirect 

^Skyrme "^Coulomb t 

CPPP 
Skyrme ' 



--DI ' 



(A2a) 
(A2b) 
(A2c) 



making explicit the power of the density matrices enter- 
ing a given term. Let us now specify these terms more 
explicitly. 
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1. The Skyrme energy functional 

We restrict ourselves here to those terms of the Skyrme 
EDF depending on time-even densities and currents that 
contribute to the ground states of even-even nuclei in SR 
and MR-PNR calculations. Also, the functional given 
below corresponds to the particular Skyrme interaction 
SIII used throughout this paper. For SIII, there arc 
no density-dependent coupling constants, but the energy 



functional can be divided into a bilinear ^g^y^jj^^ and a 
trilinear term Sg^y^.^^^,- The Skyrme energy functional 
is usually represented either in terms of isoscalar and 
isovector densities [82j or in terms of the total density 
and the densities of the nucleon species [11] . In the con- 
text of particle-number restoration, the most convenient 
representation separates contributions which are bilinear 
in densities of the same isospin from those that are bilin- 
ear in densities of different isospin 
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Skyrme 



CPPP 

Skyrme 



L q=p,n 

q.q' =p,n 

jd\ Y ^"'pIp^'- 



(A3) 
(A4) 



q.q' =p.i 

qi^q' 



The A^^ and B^^ denote the coupling constants, ^"^ none 
of which is density dependent for SIII. In the canonical 
basis, the local densities entering the energy functional 

mm 

are given by 

r,(r) = 2^ [y^lirq)] ■ [V<^.(rg)] p^, (A5) 

and denote, for the isospin q ~ n,p, the matter density, 
the kinetic density and the spin-orbit current, respec- 
tively. The operator a is the vector built out of the three 
cartesian Pauli matrices. The density matrix is ei- 
ther given by Eq. (O for the SR EDF, or by Eq. USD or 
([38)1 for the PNR MR EDF. One can see from the ex- 
pressions given above that any local density fq{r) can be 
written as: 

/,(r)^2^T^/^(rq)p^^, (A6) 

where / G {p, r, J} and where the explicit form of each 
W/f (rq) can be easily extracted from Eq. (jASp . This will 



Superscripts //' and ///' used on the r.h.s. of Eqs. ||A3IA4| | 
refer to the local densities that appear in the functional, while 
the superscripts pp, kk, ppp, ... on the l.h.s. of Eqs. ||A3| |. ||A4| |. 
and 1 IA8I 1 correspond to the powers in the density matrices. 



facilitate the construction of the matrix elements needed 
to evaluate the correction E^r- 



2. The Coulomb energy functional 

The standard Coulomb energy functional that is used 
in connection with most parameterizations of the Skyrme 
energy functional is given by 

The proton density entering Eq. (jA7p is calculated as 
described in the preceding section. The energy func- 
tional (|A7|) provides the textbook example of an energy 
functional that is not self-interaction free [2^ . 

The Coulomb exchange term in the Slater approxi- 
mation, represented by the second term on the r.h.s. of 
Eq. (|A7|) , resembles the density-dependent terms of mod- 
ern parameterizations of the Skyrme functional. As, at 
present, we do not have a correction scheme for terms 
depending on non-integer powers of the density, we drop 
it and consider the direct term only in the present work. 
Concerning absolute binding energies, the Coulomb ex- 
change term is the smallest of all contributions to the 
energy functional for nuclei and states considered here; 
it does not exceed 2 % of the total binding energy even 
in very heavy nuclei with a strong Coulomb field. What 
is even more important for the present study is that its 
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value changes also by at most 2 % when deforming a nu- 
cleus; its influence on potential energy surfaces is smaller 
than what can be resolved in the plots shown in Sec. IVII 



3. Pairing energy functional 

For pairing, we choose a local energy functional de- 
duced from a simple delta interaction (DI) , often referred 
to as "volume pairing" 



-DI 



(A8) 



More elaborate parameterizations of the pairing energy 
functional are frequently used in the literature. When 
enforcing time-reversal invariance as done here, the local 
pair densities entering the pairing functional are related 
to the pairing tensor through 



Pi 



(r) ^ 2^<^(r<j), 



M>0 



(A9) 



(AlO) 



At>0 



"MM 



and Kfi^i^^ * are given by Eqns. ([6]) and ^ for 
SR-EDF calculations with (f' = ip, and by Eqns. pT]) 
and (HH), or Eqns. ([Ml) and (gO]), respectively, for PNR 
MR-EDF calculations. In the case of SR EDF and PNR 
MR EDF calculations, W^^(rg) and Wl^^{rq) are equal 
and given by 

o-=±l 

(All) 

and represent the spin-singlet part of the two-body wave 
function. This does not hold for other MR EDF calcula- 
tions. The notation a = ±1 denotes the spinor compo- 
nent with spin projection ±1/2. The functions Wl^^{rq) 
and W^p(r(7) incorporate a cutoff to regularize the 
pairing problem, which is otherwise divergent in a vari- 
ational calculation. In the SR calculations, we use the 
smooth phenomenological cutoff proposed in Ref. [sl ]. 
while in the PAV-PNR MR calculations it is set to = 1. 



APPENDIX B: CORRECTION TERM 
1. Bilinear terms 

a. Matrix elements 

We focus here on the case where the system is time- 
reversal invariant, which leads to 



"mm 



W-- 



(Bl) 



for the time-even densities contributing to the Skyrme 
and Coulomb functionals. There is a minus sign in the 



l.h.s. of Eq. (jBip when considering time-odd ones that 
we do not have to take into account here as the corre- 
sponding contributions from the two states (/i, fl) cancel 
out both in the total energy and in the correction given 
by Eq. ((29)) . For the state-dependent function entering 
the pair density we have 

i^l) - <m (r?) = -.9m W^m'm (r?) ■ (B2) 
For the SIII energy functional, the matrix elements that 
match the definition of the bilinear part as given by 
Eq. ([9]) read as 



^2 jd^Y. (1-9) W£irq) , (B3) 

fj' 

where the sum over /, /' runs over all terms appearing 
in Eq. (|A3[) . The quasi-local form of the Skyrme energy 
functional simplifies the construction of the matrix ele- 
ments v'^l^i, in two ways: on the one hand, they involve 
one triple integral only, and on the other they contain 
products that are separable in p and v. This is of great 
help from the numerical point of view when coding the 
correction to the PNR energy as defined by Eq. (|29p. 

The situation is different for the direct Coulomb term. 
Indeed, the corresponding matrix elements (not antisym- 
metric as Coulomb exchange was dropped all together) 
are not separable 



2e' 



^3 3 ,WP^{vp)WSA^'p) 



(B4) 



and they involve a six-fold integral. This considerably 
complicates their calculation compared to the matrix el- 
ements of the Skyrme functional. Instead, the Poisson 
equation for the Coulomb potential generated by the 
source W^^{rp) 



{v) ^ -Aire ^W'^^iyp) , 



(B5) 

is solved first using boundary conditions constructed 
from the lowest-order terms in a multipole expansion of 
the state-dependent field W^^^ (rp) , and then the Coulomb 
energy of the other density in this field is obtained as 



2e' d'rWPAvp)U,A^). 



(B6) 



For all but very light nuclei, the calculation of the correc- 
tion is much more costly than the calculation of the PNR 
direct Coulomb energy itself, as the correction J7^^(r) 
has to be determined for each single-particle state solv- 
ing Eq. (jBSp . while for the total Coulomb energy the 
Coulomb potential has to be determined only for the 
summed up total charge density. However, J7^^(r) en- 
tering the correction is independent of the gauge an- 
gle, while the Coulomb potential has to be determined 
for each gauge angle when calculating the total PNR 
Coulomb energy. 

Last, but not least, the matrix elements entering the 
pairing functional are given by 



v'-'^ - = 4 



d'rAPPWl:;{rq)wUrq) 



(B7) 
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Correction 



which we introduce the shorthand notation 



Let us now write the spurious contribution £^G that 
must be removed from the MR-PNR energy, defined by 
Eq. for the functional introduced in Appendices lA 11 
IXtlandfXl 

The spurious contributions only originate from inter- 
actions between particles of the same isospin. All contri- 
butions from the bi-linear part of the energy functional 
to the correction contain the same occupation factor, for 



g-.»vJV ^g2iv _ ly 



2tt6 



(B8) 



v>0 



Hence, we obtain 



J 



cN 



CG 



Ai>0 



E 



(B9) 



where it is understood that the Coulomb term only contributes to the sum over proton pairs. In the MR-PNR code, 
the calculation of Eq. (jB9p constitutes an effort similar to the evaluation of a local one-body operator, as it can be 
reduced to a single sum over half of the single-particle states adding up a local function in space that is integrated 
over afterwards. 



Trilinear terms 



a. General expression 



We have restricted ourselves here to the special case of the Skyrme SIII functional. The zero-range three-body force 
that it originates from has the particular property that it gives an energy functional composed of terms which are 
bilinear in densities of one isospin times a density of the other isospin. The absence of terms trilinear in densities of 
one isospin greatly simplifies the correction term (see Paper I), which reduces to 
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where {N\ = N, A^^ = Z) or (A^^ = Z, = N) depending on the isospin of the states {iJ,,fJ,). 



Matrix elements 



(BIO) 



The matrix elements of the trilinear term appearing in the SIII Skyrme functional are given by 

v;Z^.X = 6 / d'rA"'"' W^^irq,) W^^ivq^) W^Avq^.) . 



(Bll) 



Correction 



Finally, the spurious term to be removed from the trilinear part of the SIII Skyrme functional is 
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where the last term in square brackets [■ • • ] is nothing but the particle-number projected local density of nucleons 
with isospin ^ g^;. 



M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. 
Phys. 75, 121 (2003). [31 
P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 
(1964). [32 
R. M. Dreizler and E. K. U. Gross, Density Functional 
Theory: An Approach to the Quantum Many-Body Prob- [33 
lem (Springer- Verlag, Berlin, 1990). [34^ 
R. G. Parr and W. Yang, Density-Functional Theory of [35 
Atoms and Molecules (Clarendon, Oxford, 1989). [36 

C. Fiolhais, F. Nogueira, and M. Marques, eds., A Primer 

in Density Functional Theory, vol. 620 of Lecture Notes [37 

in Physics (Springer, Berlin and Heidelberg, 2003). 

W. Koch and M. C. Holthausen, A Chemist's Guide [38 

to Density Functional Theory (Wiley- VCH, Weinheim, [39 

2001). 

A. Nagy, Phys. Rep. 298, 1 (1998). [40 
W. Kohn, Rev. Mod. Phys. 71, 1253 (1998). 

J. Engel, Phys. Rev. C 75, 014306 (2007). 

B. Giraud, B. K. Jennings, and B. R. Barrett, Phys. Rev. [41 
A 78, 032507 (2008). 

B. Giraud, Phys. Rev. C 77, 014311 (2008). [42 
P.-G. Reinhard and W. Otten, Nucl. Phys. A 420, 173 [43 
(1984). 

S. Aberg, H. Flocard, and W. Nazarewicz, Annu. Rev. [44 
Nucl. Part. Sci. 40, 439 (1990). 
W. Nazarewicz, Prog. Part. Nucl. Phys. 28, 307 (1992). [45 
W. Nazarewicz, Int. J. Mod. Phys. E2, 51 (1993). 
P. G. de Gennes, Superconductivity of Metals and Alloys [46 
(Benjamin, New York, 1966). 

P.-H. Heenen, P. Bonche, J. Dobaczewski, and H. Flo- [47 
card, Nucl. Phys. A561, 367 (1993). 
M. Anguiano, J. L. Egido, and L. M. Robledo, Nucl. [48 
Phys. A696, 467 (2001). 

M. Samyn, S. Goriely, M. Bender, and J. M. Pearson, [49 
Phys. Rev. C 70, 044309 (2004). [50 
T. Niksic, D. Vretenar, and P. Ring, Phys. Rev. C 74, 
064309 (2006). [51 
L. M. Robledo, Int. J. Mod. Phys. E16, 337 (2007). 
N. Tajima, H. Flocard, P. Bonche, J. Dobaczewski, and [52 
P.-H. Heenen, Nucl. Phys. A542, 355 (1992). 
F. Donau, Phys. Rev. C 58, 872 (1998). 

D. Almehed, S. Frauendorf, and F. Donau, Phys. Rev. [53 
C63, 044311 (2001). 

J. Dobaczewski, W. Nazarewicz, P. G. Reinhard, and [54 
M. V. Stoitsov, Phys. Rev. C 76, 054315 (2007). [55 
T. R. Rodriguez and J. L. Egido, Phys. Rev. Lett. 99, 
062501 (2007). [56 
T. Duguet, M. Bender, K. Bennaceur, D. Lacroix, 
and T. Lesinski (2008), submitted to Phys. Rev. C, 
arXiv;0809.2049. 

D. Lacroix, T. Duguet, and M. Bender (2008), submitted [57 
to Phys. Rev. C, arXiv:0809.2041. [58^ 
J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 
(1981). [59 
S. Stringari and D. M. Brink, Nucl. Phys. A304, 307 [60 



(1978). 

P. Bonche, J. Dobaczewski, H. Flocard, P.-H. Heenen, 
and J. Meyer, Nucl. Phys. A 510, 466 (1990). 
P.-G. Reinhard and C. Toepffer, Int. J. Mod. Phys. E3, 
435 (1994). 

N. Onishi and S. Yoshida, Nucl. Phys. 80, 367 (1966). 
R. Balian and E. Brezin, Nuovo Cimento 64, 37 (1969). 

G. C. Wick, Phys. Rev. 80, 268 (1950). 

P. Ring and P. Schuck, The Nuclear Many-Body Problem 
(Springer- Verlag, New- York, 1980). 

J. Blaizot and G. Ripka, Quantum Theory of Finite Sys- 
tems (MIT Press, Cambridge, Massachusetts, 1986). 

H. J. Mang, Phys. Rep. 18, 327 (1975). 

P. Ring, R. Beck, and H. J. Mang, Z. Phys 231, 10 
(1970). 

A. L. Goodman, in Advances in Nuclear Physics, edited 
by E. V. J. W. Negele (Plenum Press, New York, 1979), 
vol. 11, p. 263. 

D. Vautherin and D. M. Brink, Phys. Rev. C 5, 626 
(1972). 

J. Decharge and D. Gogny, Phys. Rev. C 21, 1568 (1980). 
P.-G. Reinhard and H. Flocard, Nucl. Phys. A 584, 467 
(1995). 

E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and 
R. Schaeffer, Nucl. Phys. A635, 231 (1998). 

M. Bender, J. Dobaczewski, J. Engel, and 

W. Nazarewicz, Phys. Rev. C 65, 054322 (2002). 

S. A. Fayans, S. V. Tolokonnikov, E. L. Trykov, and 

D. Zawischa, Nucl. Phys. A676, 49 (2000). 

M. Anguiano, J. Egido, and L. Robledo, Nucl. Phys. 

A683, 227 (2001). 

D. Vretenar, A. V. Afanasjev, G. A. Lalazissis, and 
P. Ring, Phys. Rep. 409, 101 (2005). 
J. A. Sheikh and P. Ring, Nucl. Phys. A665, 71 (2000). 
M. Anguiano, J. L. Egido, and L. M. Robledo, Phys. 
Lett. B545, 62 (2002). 

J. A. Sheikh, P. Ring, E. Lopes, and R. Rossignoli, Phys. 
Rev. C 66, 044318 (2002). 

M. V. Stoitsov, J. Dobaczewski, R. Kirchner, 
W. Nazarewicz, and J. Terasaki, Phys. Rev. C 76, 
014308 (2007). 

T. R. Rodriguez, J. L. Egido, and L. M. Robledo, Phys. 
Rev. C 72, 064303 (2005). 

G. Ripka and R. Padjen, Nucl. Phys. A132, 489 (1969). 
C. D. Siegal and R. A. Sorensen, Nucl. Phys. A184, 81 
(1972). 

M. Bender and T. Duguet, Int. J. Mod. Phys. E16, 222 
(2007), The correction used is this paper is not exactly 
the same as the one used here and contains a part of the 
true multi-reference self-interaction and self-pairing. 

B. F. Bayman, Nucl. Phys. 15, 33 (1960). 

B. Banerjee, H. J. Mang, and P. Ring, Z. Phys A215, 
366 (1973). 

T. Duguet (2006), unpublished. 

C. A. Ullrich, P.-G. Reinhard, and E. Suraud, Phys. Rev. 



37 



A 62, 053202 (2000). 

[61] C. Legrand, E. Suraud, and P.-G. Reinhard, J. Phys. B 
35, 1115 (2002). 

[62] A. Ruzsinsky, J. P. Perdew, G. I. Csonka, O. A. Vydrov, 
and G. E. Scuseria, J. Phys. Chem. 126, 104102 (2007). 

[63] E. Engel, Lecture Notes in Physics 620, A primer in Den- 
sity Functional Theory (Springer, 2003). 

[64] S. Kiimmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008). 

[65] K. Rutz, M. Bender, P.-G. Reinhard, and J. A. Maruhn, 
Phys. Lett. B468, 1 (1999). 

[66] M. Bender and P.-G. Reinhard (1999), unpubhshed dis- 
cussion notes. 

[67] M. Bender, H. Flocard, and P.-H. Heenen, Phys. Rev. C 

68, 044321 (2003). 
[68] T. Duguet, M. Bender, P. Bonche, and P.-H. Heenen, 

Phys. Lett. B559, 201 (2004). 
[69] M. Bender, P. Bonche, T. Duguet, and P.-H. Heenen, 

Phys. Rev. C 69, 064303 (2004). 
[70] M. Bender, G.-F. Bertsch, and P.-H. Heenen, Phys. Rev. 

C 73, 034322 (2006). 
[71] M. Bender and P.-H. Heenen, Nucl. Phys. A713, 390 

(2003). 

[72] A. P. Severyukhin, M. Bender, H. Flocard, and P.-H. 
Heenen, Phys. Rev. C 75, 064303 (2007). 



[73] M. Bender, P. Bonche, and P.-H. Heenen, Phys. Rev. C 

74, 024312 (2006). 
[74] K. Dietrich, H. J. Mang, and J. H. Pradal, Phys. Rev. 

135, B22 (1964). 
[75] V. N. Fomenko, J. Phys. (G.B) A 3, 8 (1970). 
[76] T. Biirvenich, D. G. Madland, J. A. Maruhn, and P.-G. 

Reinhard, Phys. Rev. C 65, 044308 (2002). 
[77] M. Beiner, H. Flocard, N. V. Giai, and P. Quentin, Nucl. 

Phys. A238, 29 (1975). 
[78] M. Baldo, P. Schuck, and X. Vifias, Phys. Lett. B663, 

390 (2008). 

[79] H. Flocard and N. Onishi, Ann. Phys. (NY) 254, 275 
(1997). 

[80] H. Zdunczuk, J. Dobaczewski, and W. Satula, Int. J. 

Mod. Phys. E 16, 377 (2007). 
[81] M. V. Stoitsov et al. (2007), unpublished. 
[82] E. Perlinska, S. G. Rohozinski, J. Dobaczewski, and 

W. Nazarewicz, Phys. Rev. C 69, 014316 (2004). 
[83] P. Bonche, H. Flocard, and P.-H. Heenen, Nucl. Phys. 

A467, 115 (1987). 
[84] P. Bonche, H. Flocard, P.-H. Heenen, S. J. Krieger, and 

M. S. Weiss, Nucl. Phys. A443, 39 (1985). 



